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ABSTRACT 

Analyses  of  electric  power  systems  ordinarily  include  only  a 
small  portion  of  the  entire  network.  Effects  external  to  the  network 
of  interest  are  modeled  as  leads  or  generators.   A  method  is  proposed 
for  identifying  an  equivalent  model  of  the  external  network  without 
taking  any  measurements  in  the  external  network  itself.   The  method  is 
simple  and  could  easily  be  implemented  to  supplement  state  estimation 
techniques. 

A  linear  load  flow  model  which  relates  real  power  to  voltage  phase 
angle  is  developed  for  analyzing  the  entire  network  and  from  this  two 
models  are  developed  for  relating  measurements  within  the  system  of 
interest  to  the  parameters  of  the  equivalent  model.   It  is  found  that 
there  are  two  components  to  the  equivalent  model.   One  is  a  fixed  struc- 
ture and  the  other  is  an  equivalent  power.   Since  both  models  for 
identifying  the  equivalent  network  are  in  the  foi*m  of  an  input,  output 
relation  with  an  unknown  additive  disturbance,  a  general  solution  for 
identifying  linear  static  systems  is  found  subject  to  certain  conditions. 
Applying  the  solution  to  the  two  models  for  identifying  the  equivalent 
network,  it  is  found  that  one  of  the  models  has  an  inherent  tendency  to 
produce  consistently  biased  estimates  of  the  equivalent  model  parameters. 
However,  that  same  model  has  the  advantage  that  it  uses  data  for  its 
input,  output  that  is  readily  available  and  relatively  accurate. 

A  simulation  made  to  test  the  proposed  method  is  discussed  for  a 
sample  system  and  the  results  agree  well  in  both  the  actual  parameters 
identified  and  predicted  error.  Sevpral  important  points  are  discussed 
with  a  view  toward  implementation. 
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1.   INTRODUCTION 

Analyses  of  electrical  power  transmission  networks  usually  carry 
the  implication  that  the  network  under  consideration  can  be  isolated 
from  all  other  networks.   In  actuality  this  can  seldom  be  realized 
since  the  network  is  interconnected  with  one  or  more  outside  networks, 
and  whatever  happens  to  the  network  being  analyzed  is  determined  at 
least  to  some  extent  on  how  it  interacts  with  the  others.  Therefore, 
it  would  be  convenient  if  the  effects  of  all  the  outside  networks  on 
the  network  of  interest  could  be  summarized  by  one  equivalent  network. 
In  most  cases  the  equivalent  network  would  be  considerably  smaller  than 
all  the  actual  outside  networks  that  it  would  represent.  For  instance, 
if  the  New  England  network  were  interconnected  with  the  rest  of  the 
United  States  at  a  half  dozen  points,  as  far  as  the  New  England  network 
would  be  concerned,  the  entire  rest  of  the  United  States  could  be 
shrunk  to  an  equivalent  network  having  a  half  dozen  buses.  A  method  is 
proposed  for  finding  the  equivalent  network  which  could  be  used  for 
contingency  planning  and  state  estimation  techniques  L.3  J,  I  4  J,  L5  Jj 
f8],r9],  and  [12]  . 


-8- 

1.1  General  Problem:   Identification  of  Linear,  Static  Systems 

When  analyzing  a  system,  the  usual  approach  is  to  surround  the 
system  of  interest  with  a  conceptual  boundary.  That  within  the  boundary 
is  idealized  in  some  respect  to  emphasize  certain  features.  Further, 
the  idealized  system  is  allowed  to  communicate  with  the  outside  only 
through  certain  conceptual  ports.  This  idealized  representation  of 
actuality  is  called  the  system  model  and  by  analyzing  the  model  conclu- 
sions can  be  drawn  about  the  actual  system.   If  the  contents  of  the 
idealized  system  and  what  enters  through  the  ports  are  completely  known, 
then  what  happens  to  the  idealized  system  can  be  found  at  least  in 
principle.   If  the  model  is  static,  then  only  vhat  currently  enters  is 
necessary  to  know  the  current  state  of  the  system.   In  such  a  case,  if 
only  knowledge  of  the  current  state  is  required,  then  the  effect  of  the 
system  outside  the  boundary  is  summarized  by  what  enters  through  the 
ports.  However,  to  know  what  will  happen  in  the  future,  it  is  necessary 
to  know  what  will  enter  in  the  future,  and  it  may  be  that  what  enters 
through  the  ports  is  at  least  partly  determined  by  how  the  system  inter- 
acts with  the  outside.   If  everything  outside  the  system  is  also  ideal- 
ized, then  everything  of  interest  can  be  modeled  by  the  two  idealized  sys- 
tems. The  first,  the  original  system  model,  represents  what  is  within 
the  boundary.  The  second,  the  external  system  model,  represents  all 
that  is  not  represented  by  the  first.   Further,  the  two  models  can  only 
interact  through  the  originally  specified  ports.  The  external  system 
model  may  be  rather  extensive  and  it  might  be  desirable  to  find  an 
equivalent  external  model  which  would  interact  with  the  original  system 
model  in  the  same  way. 
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The  process  of  finding  the  equivalent  system  model,  called  identi- 
fication here,  is  investigated  for  a  class  of  linear,  static  system 
models.   More  specifically,  it  is  assumed  that  the  model  consists  of  a 
structure  which  remains  constant  and  two  types  of  variables,  inputs  (or 
causes),  and  outputs  (or  effects).  The  mathematical  notation  using 
matrices  is  simply  y  =  C  u  where  u  represents  an  input  vector,  y  repre- 
sents an  output  vector,  and  C  represents  the  structure. 

If  the  original  system  model  with  its  conceptual  ports  to  the 
outside  is  placed  in  this  framework,  then  it  might  be  possible  to 
separate  what  exists  at  or  passes  through  the  ports  into  the  same  two 
types  of  variables.  For  the  model  of  a  physical  system,  the  product  of 
the  two  variables  frequently  represents  a  generalized  power  passing 
through  the  ports.   An  example  would  be  the  voltage  level  existing  at 
an  electrical  terminal  and  the  current  passing  through  the  terminal. 
By  combining  inputs  from  the  ports  with  inputs  internal  to  the  system, 
the  outputs  can  be  evaluated.   In  the  same  way  knowledge  of  inputs  at 
the  ports,  inputs  within  the  external  system  model,  and  the  structure  of 
the  external  model  allow  outputs  of  the  external  model  to  be  found.  That 
is,  each  model  may  be  analyzed  independently  of  the  other  by  knowing 
what  passes  through  the  ports.  However,  in  order  to  predict  future 
outputs  of  the  original  system,  knowledge  of  how  the  two  models  interact 
is  necessary.  That  is  the  purpose  of  identifying  the  equivalent  external 
system  model. 

The  inputs  and  outputs  can  be  divided  into  those  belonging  to  the 
original  system,  those  of  the  external  system,  and  those  common  to  both. 
Two  general  approaches  to  identifying  the  equivalent  model  are  investi- 
gated.  One  involves  combining  both  the  original  and  external  models 
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and  then  identifying  the  new  model  using  all  the  inputs  not  in  the 

external  model  and  only  the  outputs  common  to  both.  Another  approach 
is  to  use  the  inputs  and  outputs  at  the  ports  or  terminals  to  identify 
the  equivalent  external  system  separately.   Both  approaches  are  analo- 
gous to  a  black  box  which  contains  both  the  external  model  structure 
and  inputs  belonging  to  the  external  model.   Inputs  from  the  original 
system  enter  the  black  box  and  outputs  are  those  which  are  common  to 
both. 
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1.2     Specific   Problem:     Electrical  Power  Transmission  Networks 

Ordinarily  when  analyzing  electrical  power  transmission  systems,   an 
arbitrary  boundary  is  drawn  about  some  part  of  the  entire   network  and 
that  part  of  the  system   is  studied   assuming  power   flowing   in  through 
the  transmission  lines   on  the  boundary  is   known.     The  specific  problem 
investigated   is  that  of  determining  how  an  interconnected  electrical 
power  system  affects    a  particular  subsystem.     The  subsystem  will  be 
named  "Own  System*   and  referred  to  as  OS   for  convenience.     The  rest  will 
be   named   "External  System"   and  referred  to  as  XS.     Together  03   and  XS 
form  the  entire  system  named    "Whole  System"   or  V/S .     The  External  System 
may  be  the  entire  outside  world  as  viewed  from  one  region.      Or  OS  may 
be   a  higher  voltage  transmission  network  looking  down  into  a  more  com- 
plex lower  voltage  distribution  network.     XS  may  also  be  a  relatively 
extensive  region  where  power  and  voltage  measurements   are   not  made  which 
is  surrounded  by  OS  where  measurements   are  made.     To  aid  in  analysis 
notation  OS   is   further  subdivided   into   "Internal  System"   or    IS  which 
has    no   immediate   connections   to  XS   and   "Boundary  System"  which  has 
immediate   connections   to  XS.     What   is   to  be  determined    is    an  equivalent 
model  for  XS   named   "Equivalent  System"  or  ES  which  will  affect  03  the 
s  ame  as  XS . 

It  will  be   assumed   that  measurements    are   available   in  03  but   not   from 
XS.      There   are  two   aspects   which  set  this   particular   identification  prob- 
lem apart   from  many  others.      One   is  that  the   identification  must  be 
primarily  passive.      It  will   not  be  possible   to  make   any  great  manipula- 
tion of  input  signals    (power   flows)   that  unduly  upset  the  system  or  do 
not   conform  to  power  consumption  requirements.      Also   it   can  be   expected 
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that  power  requirements  in  XS  will  be  about  the  same  as  in  OS,  and 
since  the  unknown  inputs  in  XS  will  cause  unknown  changes  in  the 
outputs  of  OS,  "signal  to  noise  ratios"  will  be  about  one  to  one. 


-14- 

1>3 Contents  of  Thesis 

The  basic  conventional  model  which  is  used  for  analyzing  networks 
will  be  developed  in  Section  2.  Variables  of  interest  for  this  model 
are  real  power,  reactive  power,  voltage  magnitude,  and  the  relative 
voltage  phase  angle.  The  model  is  nonlinear  and  involves  cross  coupling 
of  all  variables.  However  it  will  be  shown  that  under  certain  conditions 
the  effect  of  voltage  magnitude  changes  on  real  power  is  much  less  than 
the  effect  due  to  phase  angle  changes,  so  that  an  approximate  linear 
model  can  be  developed  which  relates  real  power  to  voltage  phase  angle. 
The  linear  model  is  not  conventional  but  it  is  not  new  either.   One  of 
its  chief  advantages  is  a  very  significant  reduction  of  computation 
requirements.  But  even  though  a  linear  model  will  be  used  for  identi- 
fying the  model  for  Equivalent  System,  it  is  still  possible  to  incor- 
porate the  identified  equivalent  linear  model  in  the  nonlinear  model  far 
Cwn  System.   In  Section  3  the  linear  model  for  External  System  is 
reduced  to  an  equivalent  model.   Then  two  models  for  use  in  identifying 
Equivalent  System  are  developed  using  the  equivalent  model.   One  of  the 
models  is  attractive  from  the  point  of  view  that  it  only  requires  obser- 
vations at  the  tie  lines  joining  OS  to  XS,  but  it  x^ill  turn  out  that 
there  are  inherent  problems  with  this  model  due  to  the  correlation 
among  variables  at  the  tie  lines.  Both  models  will  be  placed  in  a 
general  form  and  in  Section  4.  the  solution  to  the  general  problem  will 
be  found  subject  to  certain  conditions.   The  solution,  which  has  a  very 
simple  form,  is  derived  from  both  an  assumed  mathematical  probability 
model  for  the  disturbances  and  the  method  of  weighted  least  squares. 
Error  analysis  equations  are  also  derived  to  study  the  source  of  errors 
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and  find  a  measure  of  confidence  in  the  identified  parameters. 

Applying  the  solution  of  the  general  problem  to  the  two  specific 
identification  models  in  Section  5»  the  effects  of  the  conditions 
placed  on  the  solution  are  analyzed.   One  shortcoming  of  the  more 
easily  implemented  model  is  discussed  in  terms  of  independent  inputs. 
A  simulation  made  to  test  the  proposed  methods  of  identification  is 
discussed  in  Section  6,  and  the  results  show  excellent  correspondence 
with  the  theory.  Errors  in  the  identified  parameters  are  within  pre- 
dicted accuracy.  Finally,  in  Section  7  several  aspects  of  the  specific 
problem  are  discussed  which  involve  increasing  identification  accuracy 
and  verification  of  the  identified  model.   It  is  also  pointed  out  how 
models  for  identifying  the  equivalent  system  can  be  valuable  by  them- 
selves for  the  purpose  of  network  model  reduction. 
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2 .  E IECTR IC  P.  L  POWER  TRA  IS  K  JSS 10  N  MODE  IS 

2gl Ibnlinear  load  Flow  Model 

In  reality  an  electrical  power  system  consists  of  a  large  variety 
of  devices,  many  capable  of  storing  energy  with  time  constants  and 
natural  frequencies  spanning  a  very  wide  range.  However  for  the  pur- 
pose of  computing  power  flows  in  sinusoidal  steady  state  only  a  few 
types  of  components  need  be  considered.  The  most  important  of  these  are 
the  transmission  lines  which  transmit  power  frcm  points  of  generation 
to  loads,  the  transformers,  and  the  buses  which  form  connecting  points 
or  nodes  for  transmission  lines,  transformers,  generating  stations,  and 
loads.  The  values  of  interest  are  the  power  flows  and  bus  voltages. 
Although  the  transformers,  transmission  lines,  loads,  and  sources  are 
capable  of  storing  energy,  when  the  system  is  in  sinusoidal  steady  state, 
the  usual  approach  is  to  treat  energy  storage  devices  as  complex 
impedances  or  admittances  and  the  sources  and  loads  as  injections  of 
complex  power.  Kence  although  the  system  is  oscillating,  it  can  be 
analyzed  as  a  static  network. 

A  transmission  line  is  of  course  a  distributed  parameter  device. 
However  as  with  most  such  devices,  it  is  possible  to  model  it  as  a  set 
of  lumped  elements  provided  the  physical  size  is  small  compared  to  the 
wavelength  of  the  power  transmitted.  Since  the  wavelength  of  60  hertz 
alternating  current  is  roughly  the  distance  between  New  York  and  San 
Francisco,  while  most  transmission  lines  are  far  less,  lumped  element 
models  are  quite  valid  for  this  purpose.  The  pi  model,  common  for  this 
type  of  distributed  parameter  device,  is  used  here  and  consists  of  a 
resistance  and  Inductance  in  series  with  shunt  capacitance  at  each  end. 
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Figure  2.1   Pi  Transmission  Line  Model 
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Figure    2.2      Electrical   Power  Network   Lumped   Model 
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Since   the  system  studied   is   assumed  to  be   in  sinusoidal  steady 

state,   the  capacitances    and   inductance  can  be   represented    as   admit- 
tances.     A  bar  above   a  symbol  will  denote   a  complex  quantity    (e.g. 

jDi 
E.=  E     e  =  E     cos  D.   +   jE     sin  D.).      An  underlined  symbol  will 

denote   a  matrix  or  vector.      Let 

-J'0ik  1 

7ik   "  yik  e  "     R.k  +   jcoL.k 

be  the  complex  admittance  of  the  series  resistance  and  inductance,  and 
let 

ysik  =  *°  — 

be  the   admittance  of  the  shunt  capacitance.     The  current  entering  the 
transmission  line   at  the   i       node   is      IL.    =  E.  ys,,    +    (E.    -  E,  )  y.,    . 

X  1  IK  x  iC  IK 

— *    — *  — 

Complex  power  is   P  +   jQ  =  E    I     (I     is   the  complex  conjugate  of  I;  so 
the  power  entering  the  transmission  line   from  the   i       node   is 

The   network  model  consists   of  admittances  which  represent  the 
transmission  lines   and   transformers   connected  to  buses    or  nodes.      Also 
connected  to  each   node   is    a  current  source  which  represents   the   load 
or  generator   at  the  bus.      Power   injections    into  the  buses    are  positive 
for  generators   and   negative  for  loads.      As   shown  in  Appendix  A,    a  bus 
admittance  matrix  can  be   formed  which  represents   the  relation  between 
the   current  sources    at  each  bus   to  the  bus   volteges.      Let  E,         be   a 
vector  whose  elements   are   the  complex  bus   voltages,    I,         a  vector  of 
the   corresponding  values   of  the  complex  bus   current  sources,   and  Y, 


-19- 

the   coirplex  bus   admittance  matrix.      Then  the   relation  is 

---  th      -N-_ 

L    =  Y,   E,    or  for  the  i   bus,   I.  =  >]  Y..E.   where  there  are 

^bus   -bus  -bus  i  /-<      ik  k 

th 

N  buses.   Correspondingly,  the  complex  power  injected  into  the  i   bus 

by  the  current  source  at  the  bus  is 

_»  N  _  _      -jD   N     jD      -j9 
i   d  i    i  £<   ik  k    i      j^i     k      ik 

=  2E.EkY.ke    ik  *  (,.2) 

k=l 

The  load  flow  problem  is  that  of  given  bus  power  injections,  find 
the  bus  voltages  and  the  transmission  line  power  flows.   It  should  be 
noted  that  even  though  the  network  is  modeled  with  linear  elements,  the 
power  flow  equations  are  nonlinear. 
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2.2  Linear  Real  Power  -  Voltage  Angle  Model 

Under  certain  assumptions  it  is  possible  to  reduce  the  complexity 
of  the  model.  From  Equation  (2.1)  line  real  power  flow  is 


PL..    =  Real 
ik 


\jEi  yslk  +  Ei  yike  "  EiEkyike  J 

=  E .2v ik  cos  0ik  -  E iE]7ik  cos  (0lk  +  D.   -  Dk) 

PLik  "  Vik"5  "ik[Ei  "  Ek  cos(Di  "  W  EiVikSin  "ik  sin(Di  "  V 

(2.3) 

If  it  can  be  assumed  that  voltage  magnitudes  do  not  vary  much  from 
nominal  operating  points,  that  the  transmission  line  reactance  is  much 
greater  than  the  resistance  so  that  0.,  »  90°,  and  that  voltage  angle 

differences  (D.  -  B,  )  are  small,  then  when  E-  =  (E,)   .   .,  and 

l  k  '  1  1  nominal* 

E.    =   (E.)        .      _,     PL..    «    E.E  v     sin(D.    -D.  ).     For  convenience  it 
2  2  nominal '  ik  1  2  ik  i         k 

^»         «* 
will  be  assumed  that  E_    =  E     =1.0.      Using  the  approximation  for  the 

sine,  sin(D     -  D,  )    ft*  D     -  D     for  D.    -  D     small 

PLik "  yik(Di  -  V  <2-^ 

As   shown  in  Appendix  B     for  such   a  case   a  linear  load   flow  equa- 
tion,  P,         =  Y,      D        ,   relates   real  power  injections  to  voltage  angles. 
For  a  system  with    (K+l)   buses,  IV        is    a  Kxl  vector  of  bus   real  pover 

injections,  D,         is    a  Kxl  vector  of  bus  voltage   angles,    and  Y,         is   a 
J  '  -bus  -bus 

KxK  matrix  whose   elements    are    nearly  the  same   as  the  magnitude   of 

elements   of  the  comolex  matrix  Y,       .      The    (K+l)       bus    is   the  reference 

-bus 

bus  at  which  the  voltage  angle  is  specified.   Power  injection  at  the 
reference  bus  is  such  that  the  algebraic  sum  of  real  pover  injections 
of  all  buses  is  zero.   Rom  [  6  J  and  Rom  and  Schweppo  C  8  3  proposed  the 
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use   of  this    linear   load   flow  model   for  on-line  estimation  of  voltage 

angles   and  transmission  line  power  flows  while  Baughman  and  Schweppe  C2  3 

showed  how  its  results   compare  with  the   nonlinear  model. 


2.3     Discussion  of  the   Linear  Load  Flow  Model 

Although  the   linear  load   flow  model  does    not  give  the   exact  values 
of  voltage   angles   obtainable  by  the   nonlinear  model   nor  does    it   account 
for  voltage  magnitude  variations,   there   are  some  great   advantages  to 
using  it   for  real  power  flow  analyses   even  when  the  conditions   on  low 
line  resistance   and   small  voltage  magnitude  spread   are   only  mildly  met. 
Foremost   is  the  simplicity  of  calculation.     For   (K+l)  buses  the   non- 
linear model  requires   an  iterative  solution  of  a  set  of  2K  simultaneous 
nonlinear  equations.     The  linear  model  involves   K  simultaneous   linear 
equations.      If  the    Fewton-Raphson  method  is  used  to  solve  the   first, 
some  comparison  between  the  two  can  be  made.     The    tewton-Raphson 
method  for  one  iteration  is   of  the   form     Aa     =  B     fee   .     Since  calcula- 
tions  for  solutions  of  simultaneous   linear  equations   are   about  propor- 
tional to  the  square  of  the   number  of  equations,   one   iteration  for  the 
nonlinear  model  requires   about   four  times  the   calculations    for  the 
linear  model.      Also  the  B     matrix  must  be  calculated   at   each   iteration 

whereas  Y,         in  the   linear  model  remains    fixed   for  a  given  system, 
-bus  J 

Finally,   depending  on  the   accuracy  required,   the    Kewton-Raphson  method 
requires    about   four  or  more    iterations.      Another   advantage   is   the  ease 
with  which  the  linear  equations   can  be  manipulated.     This   will  become 
very  important   in  transforming  the    network  for   finding   an  equivalent 
system  as    it  results    in  much   insight.      Historically,   before  digital 
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computers  were  in  wide  use,  load  flow  snalyses  were  made  with  calculating 

boards  which  used  physical  components  to  model  the  network.  The  more 
versatile  P.-C   calculating  boards  could  account  for  voltage  magnitudes 
and  line  power  losses  but  were  expensive  and  there  were  only  about  50 
in  the  United  States  L11J.  The  linear  load  flow  model  is  somewhat 
analogous  to  the  D-C  calculating  boards  which  were  more  numerous  and 
less  expensive.  The  linear  load  flow  model  accounts  for  individual  bus 
real  power  injections  while  some  D-C  calculating  boards  used  only  one 
source. 

The  error  involved  in  using  the  linear  load  flow  model  of  course 
depends  on  how  well  the  three  conditions  are  met.  The  assumption  of 
small  voltage  angle  differences  between  the  ends  of  a  transmission  line 
usually  results  in  the  least  error.  Except  in  rather  rare  circumstances, 
angle  differences  are  within  30°  and  usually  they  are  within  10°.  Com- 
paring the  values  of  the  sine  of  an  angle  and  the  angle: 


5.ngle    (degrees 

) 

Angle    (radians) 

Sine 

10° 

fg  =  0.1745 

0.1736 

20° 

\-  =  0.349 

0.342 

30° 

f-  =  0.524 

0.500 

Even  up  to  30°  the  approximation  is  within  5%,   and  under  10°  it  is 
within  0.5*. 

The  error  involved  in  the  other  two  assumptions  can  be  analyzed 
through  the  nonlinear  equation  for  line  power  flow.   If  the  line 
conductance,  g.,  =  y,,  cos  0.,  ,  and  the  line  suseptance,  -b.,  =  -y.v.  sin  ^-?v» 
are  used  then  equation  (2.3)  is 
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PLik  ■  =1^1^!  "  Ek  "■&!  -  Vl  +  Wik  sin(Di  "  V  (2-5) 

In  order  to  make  a  rather  rough  estimate  of  errors   let   "a"  represent 
either  b.,    or  y^.     Then  when  PL.,    is   approximated  by  b.,  (D     -  D,  )   or 
y.,  (D.    -  £,)  >    for  small   angle  differences 

PLtk  *  Vik(Ei  -  V +  EiEkbik(Di  -  V 

and  the  relative  error  for  PL.,  a*   a(D.  -  D,  )  is 

ik  1         k 

i   +  •                           PL  -  PL         e 
relative  error   =     — rr =     c 

n  a(Di-V  -Eigik(Ei-  V  -  EiEkbik(Di  -  V 

Wik(Di  -  v 


a  -  EiEkbik      1  Sik  Ei  - E 


k 


D     -  D 
E.E.b..  E,      b.,      ui       ^k 

l  k  lk  k       lk 

For  E.    &  1.0  and  E      ?»  1.0,     €.    is   approximately 

a  gik     Ei   "  Ek 

6   **    b..     "  Vk         b..      D,    -  D. 

ik  lk       i         k 

If  the  voltage  drop  is  expressed  as  a  percentage,  AE  =  1C0(E.  -  E,  ) , 
and  if  angle  difference  is  expressed  in  degrees,  AB  =  57.3  (D.  -  B,  ) 

1  K 

«*  50(D.    -  D,  ) ,   the  relative  error   is   approximately 

'ik         -"-  "'  "     "ik 


(fc-vO.-iSS  *•« 


The  purpose  of  equation  (2.6)  is  to  show  how  the  greatest  portion  of  the 

error  enters  the  linear  transmission  line  power  flow  equation  (2.4-). 

The  term  r —  -  E  E,   shows  how  voltage  magnitude  changes  from  nominal 
Dik    1  K 


-24- 

1  gik  *E 

values  affect  the  error,  while  -  r —  rrr  shows  the  effect  of  resis- 

2  Q.,  hu 

tance  and  voltage  drop.  Table  2.1  compares  line  flows  as  computed 
both  by  the  nonlinear  and  linear  equations  along  with  actual  errors 
and  approximate  errors  predicted  by  equation  (2.6).   It  is  not  immedi- 
ately obvious  whether  b  ,  or  y.,  (actually  E.E  b.,  or  E.E^y  ,  where 
E  and  E,  are  nominal  values)  should  be  used  for  a  in  PL.,  »  a(D  -  D,  ). 

In  view  of  the  relatively  small  error  and  much  less  involved 
calculations  required,  and  to  preserve  insight,  the  linear  load  flow 
model  is  used  for  the  purpose  of  identifying  the  equivalent  system  model. 
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3.  MODELS  FOR  IDENTIFICATION  OF  THE  EQUIVALENT  SYSTEM 

3.1  Tie  Line  Pover  Flow  Model 

As  previously  mentioned,  tvo  approaches  to  finding  the  equivalent 
system  model  will  be  used.  The  tie  line  power  flow  model  views  the 
external  system  as  a  separate  system  and  attempts  to  find  an  equivalent 
model  by  using  the  tie  line  power  flows  and  bus  voltage  angles  where  the 
boundary  around  CS  crosses  the  tie  lines  connecting  CS  to  XS.  The  bus 
injections  and  corresponding  bus  voltage  angles  can  be  separated  into 
3  groups,  those  belonging  to  IS  (internal  System),  those  belonging  to 
BS  (Boundary  System),  and  those  belonging  to  XS  (External  System). 
Collectively  IS  and  BS  form  CS  (Ovn  System).   PI  and  DI  are  the  bus 
real  power  injection  vector  and  the  bus  voltage  angle  vector  respectively 
whose  elements  belong  to  IS.  Similarly  PB  and  DB  belong  to  BS ,  and  PX 
and  DX  belong  to  XS.  The  letters  I,  B,  and  X  will  usually  refer 
values  belonging  to  IS,  BS ,  and  XS  respectively.  Since  by  earlier 
definition  IS  and  XS  are  not  immediately  connected  to  each  other  by 
transmission  lines,  the  linear  load  flow  equations  can  be  written: 


F3 
PX 


YII       Y33        0 
YJEB'     YBB     Y3X 
0         YBX'   YXX 


HI 
DB 
DX 


(3.1) 


(The  transpose  of  a  matrix  will  be  denoted  by  YT  =  transpose  of  Y). 
The  zero  submatrices  in  the  corners  are  due  to  the  lack  of  direct 
coupling  between  IS  and  XS.  The  only  coupling  between  OS  and  XS  is 
reflected  by  the  submatrix  YBB.  Since  the  elements  of  YBB  are  values 
of  admittances  of  transmission  lines  which  have  at  least  one  of  their 
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ends   connected   to  buses    in  BS,   Y3B  can  be   separated   into  that  which 

belongs   to  OS,  YBOS,    and  that  which  belongs   to  XS,  YTL,   so  that 
YBB    =  YBOS  +  YTL.      YTL  will  be   a  diagonal  matrix  whose   elements  will  be 
the   admittances   of  the  tie   lines   connecting  OS   to  XS.      Equation   (3,1) 
can  be  written 


PB 
PX 


YTB 


YII 

YT3'     YBOS 
0 


C 

0 
0  0 


51 

DB 
DX 


0 

C~ 

51 

YTL 

YBX 

DB 

YBX' 

YXX 

DX 

or 


PI 


F3 


YII       YI3 


YTB'      YBOS 


DI 


DB 


YTL  DB  +  YBX  DX 


(3.2) 


PX   =  YBX'DB  +  YXX  DX 


Let     PTL  =  YTL  D3  +  YBX  DX 


(3.3) 
(3.  A) 


Tie   line   power   flow,   PTL,    is   a  vector   of  length  equal  to  PB   and  repre- 
sents the  power  which  flows  out  of  the  buses  belonging  to  BS   into  the 
transmission  lines   which  connect  BS  to  XS.      Solving   for  DX   in  equation 
(3.3)    sn^  substituting   into  equation    (3.-0 


DX   =   -YXX_1YBX'DB  +  YXX-1PX 


PTL  =  YTL  DB  +  YBX  (-YXX'^TBX'DB  +  yXX^PX) 
=    (YTL  -  YBX  YXX'-hfBX '  )D3  +  YBX  YXX"1PX 


Define   an  equivalent   external  system  bus    admittance  matrix,  YE OX ,    and   an 
equivalent  external  bus   power  injection  vector,   FEQX ,   such   that 
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YEQX    =  -YBX  YXX'-^BX1 


PEOX   =  AQX  PX 


AQX     =  -YBX  YXX 


-1 


(3.5) 
(3.6) 
(3.7) 


Then  equations    (3.2)   and    (3.4.)   can  be   rewritten 


PI 

_ 

P3   - 

PTL 

III 

YB" 

21 

YB' 

Y30S 

D3 

PTL  =    (YTL  +  YEQX)D3   -  PEQX 


(3.8) 


(3.9) 


Equation  (3.8)  is  simply  the  linear  load  flow  model  for  OS.  By  including 
PTL  with  PB  the  problem  can  be  solved  as  usual.  Equation  (3.9)  for  PTL 
shows  two  contributions.  The  first,  (YTL  +  YEQX)DB  is  the  part  of  PTL 
due  solely  to  the  values  of  voltage  angles  existing  at  the  buses  of  BS. 
The  second,  PEQX  is  the  equivalent  power  injections  in  XS  as  seen  by 
BS.   It  is  the  bus  power  injections  from  PX  that  are  routed  to  BS  by  the 
structure  of  XS.   It  would  be  desirable  to  know  both  YEQX  and  FEQX,  but 
even  knowledge  of  just  YEQX  would  be  valuable.   For  example,  if  the 
voltage  angles  and  line  power  flows  in  OS  were  wanted  after  a  change  in 
bus  power  injections  in  OS,  then  if  PX  (and  thus  PEQX)  did  not  change 
between  times  t,  and  t  ,  the  change  could  be  found  by  substituting 
equation  (3.9)  into  (3«8)  so  that  at  any  time  t 


Pl(t) 


FB(t)  +  PF2X(t) 


YII    YIB 


YB'   (YBB  +  YEQX) 


DI(t) 


D3(t) 


and    for  the  change  with  PEQX  (t    )    =  PECX(t   ) 


-29- 


Pl(t2)-  Pl(t1) 


P3(t2)-  ra(t2) 


II         YB 


YIB'      (YBB  +  YEQX) 


'Dl(t2)    -  Dl(t1) 


D3(t2)    -  D3(t1)J 


where  Dl(t.)    and  D3(t,)    are   found  by  measuring  PTL(t..)    and  using 
equation   (3.8).      On  the   other  hand,    if  the  voltage   angles   and   line 
power  flows   in  OS  were  wanted   after  a  change   in  the  status   of  one  or 
more   transmission  lines    in  OS,  then  if  YII,   YI3,    and  YBB   are   the   admit- 
tance matrices   after  the  change,    for   Pl(t   )    =  Pl(t,),    PB(t   )    =  PBjt..), 
and   PEQX(t   )    =  FEQX  (t^ 


YII       YIB 


Y  IB '      Y3B+YEQX 


DlCt^ 


OlCt^ 


YII       YIB 


YIB'      YBB+YEQX 


2i(t2) 


_2B(t2) 


Everything   is   known  except  Dl(t   )   and  P_3(tp)   since   again  Dl(t,)    and 
DB(t,)   can  be    found  by  measuring  FTl(t,)    and  using  equation    (3.8). 

As    formulated,    in  equation    (3.9)  E>B   is   an  input  vector   and   FTL  is 
an  output  vector  while  FEQX  might  be  thought  of  as   an  unknown  distur- 
bance.     If  a  series   of  input,   output  measurements  were   taken,   then  it 
might  be  possible  to   find  YEQX.     Equation    (3.9)   will  be   referred  to  as 
the  tie   line  power  flow  model. 
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3.2  Boundary  Bus  Impedance  Model 

The  second  approach  involves  an  impedance  matrix  instead  of  an 
admittance  matrix.   In  order  to  obtain  the  required  form,  an  equivalent 
internal  system  as  seen  by  BS  looking  into  IS   will  be  found  which  is 
similar  to  the  equivalent  external  system.   The  top  row  of  equation  (3.8) 
is 

PI  =  YII  DI  +  YJB  BB 

Solving  for  DI 

DI  =  -YII"1YJB  BB  +  YII^PI 

and  substituting  into  the  second  row  of  equation   (3.8) 

P3    -  FTL  =  YB'  (-Yir1YIB  BB  +  YIl"1?!)   +  YBCS  DB 

Using  equation  (3.9)   for  PTL  and  combining  terms 

FB  +  FZQX   -  YIB'Yir1?!  =    (YBCS  +  YTL  +  YEQX   -  Y3B  ,YH"-1TI3)DB  (3.10) 

Define    an  equivalent   internal  admittance  matrix,   YEQI,    and   an  equi- 
valent  internal  bus   power  injection  vector,   HBQI,   similar  to  YEQX   and 
FEQ2  by 

YEQI  =   -IJB'YII"1YJB  (3.11) 

FEQI  =  ApJ  PI  (3.12) 

API     =  -YJB'YII"1  (3.13) 

Then  equation   (3. 10)  can  be  rewritten. 

FEQI  +  P3  +  PEQX   =    (YEQI  +  YBB  +  YEQX)D3  (3.H) 
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PEQI  is   the  power   injection  of  PI  which   is   routed  to  the  buses   of  BS 
by  the   structure   of   IS.      YEQI  is   the  structure  of  IS   as   seen  by  BS. 
YB3  was   previously  separated   into  that  belonging  to  OS,  YBCS    and   that 
belonging  to  XS,  YTL.      If  YBCS    is    further  separated    into  that  belonging 
to    IS,  YBIS   and   that  belonging  to  BS ,   YBBS   so  that  Y3C6    =  YBIS  +  YBBS 
or  YBB   =  YBIS  +  YBBS  +  YTL.,   then  the  quantity 

(YBIS  +  YEQl)D3   -  PEQI 

is  the  power  flowing  out  of  the  buses   of  BS   into  the  transmission  lines 
connecting    IS   and  BS.      If  ES  had   no  transmission  lines   unique  to   itself, 
then  YBBS  would  be   zero.      By  rewriting  equation   (3.14.)   with 
YBB    =  YBIS  +  YBBS  +  YTL 

PEQI  +  PB  +   PEQX   =    (Y3JS  +  YSQl)PB  +  YBBS  C3  +    (YTL  +  YE0X)DB  (3-Ua) 

Equation  (3.14.a)  can  be  considered  to  be  two  matrix  Thevenin  or  Norton 
equivalent  circuits,  one  for  IS  and  one  for  XS,  coupled  to  the  circuit 
for  BS  as  illustrated  in  Figure  3.1.  If  equation  (3. 14-)  is  solved  for 
DB,   then 

DB   =    (YEQI  +  YBB  +  YEOj)'1  (F£0I  +  PB  +   PEQX) 

As   shown  in  Appendix  B,    (YEQI  +  YBB  +  YEQX)"1   is   that  part   of  the  Whole 
System    (MS)   impedance  matrix  which  relates   P3  to  DB.     Therefore  denoting 
ZBB    for  that  part   of  the   impedance  matrix 

DB   =  ZBB  (PEQI  +  PB)   +   ZBB   PEQX  (3.15) 

In  this    formulation    (PEQI  +  PB)    is   an  input  vector,  DB    Is    an  output 
vector,    and   ZBB   PEQX   is    an  unknown  disturbance  vector.      Again  if  a 
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Figure    3.1     Equivalent   I'etwork  Viewed  From  BS 
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Figure    3.2      Tie    Line    Power  Flow  I.lodel 
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Figure   3.3      Boundary  Bus   Impedance   .iodel 
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series  of  input,  output  measurements  were  taken,  it  might  be  oossible 
to  find  ZBB. 

For  comparison  purposes  and  in  order  to  be  able  to  investigate 
sources  of  approximation  error  later,  the  analogous  equations  for  the 
nonlinear  model  are  developed  in  Appendix  C. 


3.3  Separating  Own  System  Estimation  From  Equivalent  System  Identification 

The  tie  line  pover  flow  model  and  the  boundary  bus  impedance  model 
are  illustrated  in  Figures  3.2  and  3.3.   As  shown  in  Table  3.1,  both 
approaches  result  in  equations  of  the  form  z  =  H  u  +  v  where  u  is  the 
input  vector,  z  the  output  vector,  v  a  vector  of  unknown  disturbances, 
and  H  the  structure  matrix  to  be  identified.   Both  z  and  u  were  assumed 
to  be  known.   However,  since  the  state  of  OS  (i.e.  the  bus  voltages)  can 
be  estimated  by  knowing  what  power  enters  from  XS  through  the  tie  lines, 
the  problem  of  estimating  the  state  of  OS  can  be  separated  from  identi- 
fying the  parameters  of  the  equivalent  system.  VJhat  happens  to  CS 
depends  on  XS ,  but  the  effect  of  XS  is  summarized  by  the  power  flowing 
through  the  tie  lines.   Hence  existing  state  estimation  techniques  [  5  3, 
[8  3,  and  [9    ]  can  be  applied  to  obtain  z  and  u.   In  practice  there  will 
be  some  error  associated  with  metering  power  and  estimating  voltage 

angles.   In  such  a  case  let  v  be  the  error  of  outout  estimate  and  v 

~z  -  -u 

be  the  error  of  input  estimate.  Then 

(z  +  X„ )  =  H(u  +  v  )  +  v 
—   — z     —    u 

or        z=Hu+(v-v+Hv) 
-   -z u 


-34- 

Table  3.1 

Ident: 

Lfication  Models 

Genera! 

L           Boundary  Bus 

Tie   Line 

Model 

Impedance 

Power  Flow 

Model 

Model 

Output 

z 

DB 

PTL 

Input 

u 

(FB+PEQI) 

D3 

I'foise 

V 

ZBB   FEQX 

-PE^X 

Structure 

H 

Z3B 

YTL+YEQX 
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By  replacing  v  -  v  +  H  v   with  a  new  unknown  disturbance  vector  v 

z  =  H  u  +  v 

the  same  approaches  can  be  used.  The  only  difference  is  that  there  is 
generally  more  disturbance  associated  with  identifying  H. 


3.4.  Use  of  Changes  in  Variables 

Up  to  now  z,  u,  and  v  have  not  been  associated  with  time.  From 
now  on  z(t),  u(t),  and  v(t)  will  be  used  to  designate  a  set  of  z,  u, 
and  v  for  a  particular  time  t.  The  disturbance  vector  v(t)  will  be 
treated  as  a  noise  term  since  u(t)  is  analogous  to  an  input  signal  and 
z(t)  is  analogous  to  an  output  signal  corrupted  by  noise.   Ordinarily 
when  working  with  estimation  or  identification  problems,  it  is  much 
easier  if  noise  terms  have  a  zero  mean.   In  the  two  approaches  above 
neither  will  in  general  have  zero  temporal  mean  noise  vectors  v(t).   It 
is  expected  that  v(t)  will  vary  about  some  value  depending  on  the  time 
of  day,  the  season,  etc.   However,  if  the  differences  in  z(t),  u(t), 
and  v(t)  are  taken  between  two  different  times,  then  it  may  turn  out 
that  the  difference  or  change  v(n)  =  v(t  .)  -  v(t  )  will  have  a  zero 
mean.  The  general  equation  is  then 

z(n)  =  K  u(n}  +  v(n)  (3.16) 

where  z(n)  =  z^^)  -  ^^n^ 

u(n)  =  u(twl)  -  u(tn) 
v(n)  =  v(tnfl)  -  v(tn) 
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It  will  turn  out  that  there  are  other  reasons  for  using  the  difference 
in  variables  rather  than  the  actual  values  of  the  variables.  The  most 
significant  is  that  using  the  actual  values  may  result  in  subtracting 
variables  whose  magnitudes  are  very  large  compared  to  their  difference. 

From  now  on,  any  reference  to  the  previously  defined  symbols  for 
voltage  angles  or  power  will  refer  to  the  changes  of  variable  between 
two  times  rather  than  the  actual  values  of  the  variables.   Hence  the  tie 
line  power  flow  model  is 

FTL(n)  =  (YTL  +  YEQX)D3(n)  -  PEQX(n)  (3.17) 

and   the  boundary  bus    impedance  model   is 

DB(n)    =  Z3B[FEQI(iQ  +  F3(n)]    +  Z33   PEQX(n)  (3.18) 

where  PTL(n)    =  PTL(t      ,)    -  FTL(t    ) 

'       n+1         rr 

P3(n)      =  PB(t     _)  -  PB(t  ) 

—  —     n+1         —     n 

D3(n)      =  DB(t      .)"  -  DB(t    ) 

—  —     n+1         —     n 

FEQl(n)    =  FEQl(t      ,)    -  PEQl(t    ) 

P£QX(n)    =  PEQX(t      ,)    -  FEQX (t   ) 
^  *-v   n+1         *-     n 

for   N  sets   of  measurements   n  =  1,2,...,N. 
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U.      IDE  Iff  IF  2CAT  ION  OF  LINEAR,  STATIC  SYSTE>5 

4-.1  Conditions  on  the  System 

Eoth  formulations  for  identifying  the  structure  of  the  equivalent 
system  have  resulted  in  the  same  general  equation,  z(n)  =  H  u(n)  +  v(n) 
where  there  are  N  sets  of  input  vectors,  u(n),  and  output  vectors,  z(n). 
H  is  the  structure  matrix  and  v(n)  is  a  pseudo  noise  vector.  The 
general  problem  can  be  stated: 

Given:  A  set  of  input  and  output  vectors  u(n)  and  z(n),  n  =  1,2,...,N 
Find:   The  structure  matrix  H  and  statistical  information  on  v(n). 

Several  assumptions  on  u(n),  v(n),  and  H  will  be  made  and  later 

these  will  be  justified  when  the  results  are  applied  to  the  specific 

problem.   First  elements  of  H  are  assumed  to  be  independent.  Specifically 

it  is  assumed  that  H  is  not  symmetric.  This  is  certainly  not  true  for 

the  specific  problem,  but  various  justifications  will  be  given  later. 

Second,  v(n)  is  assumed  to  be  a  zero  temporal  mean  process,  or 

£/X(n)r  =  Q. .   where  E^a|  is  the  expected  value  of  a.  Third,  v(n)  is 

uncorrelated  in  time,  or  Ejv(n)  v'(m)/  =  0  for  n  /  m  (n  represents  a 

sequence  in  time).   Fourth,  the  covariance  matrix  for  v(n)  is  assumed 

to  be  of  the  form  -y~T  R,  or  Efv(n)v»  (n)7  =  R(n)  =  -i-r  R.  c(n)  can 

c(.n)  -'  v-v    '—        'J        -N  c(n;  - 

+  Vi 
be   called   a  confidence  coefficient   for  the   n       measurement  set   and 

indicates   the  relative   noise   level  for  that  set.      Larger  c(n)fs    indicate 

less   noise.      Fifth,   v(n)    is   uncorrelated  with  u(n),   or     EJv(n)u' (n)v    =  0. 

At  this   point  there   are  many  possible   approaches  to  the  problem. 

Two  will  be  used   here.      One    is    a  time   tested  mathematical   approach 

which   assumes    a  particular  probability  model   for   the   noise,    and  the 

other   is   an  engineering   approach,    also  time  tested,  which  uses    no  model 
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for  the  probability  distribution  of  the  noise.   It  will  turn  out  that 

both  result  in  the  same  solution  which  will  have  a  relatively  simple 

form.  The  two  approaches  ere: 

1.  Assume  that  the  noise,  v(n),  has  a  Gaussian  probability 

distribution  (with  a  zero  mean  and  covariar.ce  of  ,  •  <   R). 

c(n)  -' 

Then  find  the  maximum   likelihood  estimates   of  H   and  R. 

2.  Assume   no  particular  probability  model  for  the   noise   and 

A 

use  the  method  of  weighted  least  squares  to  find  H.  Then 

make  a  reasonable  estimate  for  R. 
The  fact  that  the  two  approaches  result  in  identical  solutions  is  a 
well  known  consequence  of  using  quadratic  criteria  for  optimization 
and  estimation  of  linear  systems.  However,  in  addition  to  analytical 
convenience,  it  is  also  well  known  that  in  many  instances  the  Gaussian 
probability  distribution  itself  is  a  rather  reasonable  probability 
model. 


U. 2     Maximum  Likelihood  Identification 

First  the  solution  will  be  found  for  the  maximum  likelihood  approach 

assuming  Gaussian  noise.   There  are  four  known  items:  the  set  of  input 

vectors  u(n),  n  =  1, ?,...,  N,  the  set  of  output  vectors,  z(n),  n  =  1,  ?,...,  N, 

the  corresponding  confidence  coefficients,  c(n),  n  =  1,  ?,...,  N,  and  the 

form  of  the  probability  distribution  function  for  v(n),  n  =  1,2,...,N, 

which  is  zero  mean  Gaussian,  with  a  covariance  matrix  ■  /  v  R.   Note 

c(n)  - 

that  R  itself  is  not  known.  Therefore  if  estimates  of  H  and  R  can  be 
found,  the  problem  is  solved  since  knowledge  of  R  completely  specifies 
the  statistics  of  v(n)  when  v(n)  is  Gaussian. 
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The  probability  distribution  function  for  the   Kxl  vector  v(n)    is 


r        Y       v  1  ~o       -o  vf (n)c(n)R"  v(n) 

p[v(n)]    =    [(?rr)Kc-K(n)lR\j 


2         2 

e 


The   likelihood   function  of  z(n)    for   a  given  H   and  R  when  z(n)    and  u(n) 
are   Kxl  vectors    and  H   and  R  are   KxK  matrices    is 

1 

p[z(n),u(n):H,R]    =    [(2rr)   c"K(n)\Rl]"2     e<^n)  >^n)  ^>*l 

where 

j[z(n),u(n):H,R]    =  |[z(n)    -  E  u(nj]  '   c(n)   R_1[z(n)    -  H  u(n)] 

Since  v(n)    is  uncorrelated   in  time    (Esv(n)v'  (m)j    =  0  for   n  /  m) ,  the 

joint  probability  distribution  function  for  the  set  V   =  £v(l)  ,v(2)  , . . .  ,v(N)> 

is   the  product  of  the  distribution  functions   for  v(l)  ,v(2). .  .  ,v(N). 

p(V)    =  p[v(l)]    p[v(2)]    ...   pfe(N)] 

N  A         42  v'Cn^Cn^vCn) 

-    [(^)NK|R|ME    c-K(n)]    'xe2-1 
n=l 

Let  p(Z:H,R)  be  the  likelihood  function  of  the  sets 

Z  =  f  z(l),z(2),...,z(N)J,  and  U  =fu(l),u(2),...,u(N)l  for  a  given  H 


and  R.   Then 


(Z,U:H,R)  =  fa)™   lR|N  H  c"K(n)]    e'^^S)  U-U 

n=l 


where 


J(Z,U*H,R)  =  \    Zj>(n)  "  H  y(n^l  c<n)  5  [z(n)  "  2  u(n)] 
"  n=l 

A.         A 

Let  H  and  R  be  the  values  of  H  and  R  which  maximize  the  likelihood 
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function  when  Z   and   U  are  the   actual  measured   sets.      As   developed   in 
Appendix  D  these  maximum,  likelihood  estimates   are 


A 

H  = 


N  ITn  1-1 

2    c(n)z(n)u'(n)       >   c  (n)u(n)u»  (n)  U.2) 

n=l  J 1  n=l  J 


and 

N 

I 

provided  the  vectors  u(n),  n  =  1,2,... 9N  are  such  that  the  inverse  of 


N 
R=|2   c(n)[z(n)   -  H  u(n)][z(n)   -  H  u(n)]  '  (4.3) 


r  n  -1 

1  ^.    c (n)u(n)u' (n)l    exists. 
ln=l  J 


Two  things   should  be   noted.      First,    it  is    not   necessary  to  know 

A  A 

R   in  order  to  find  H.     This   is  significant   in  that  it  is   not   necessary 

A  A 

to  solve   for  R  and  H  simultaneously  which  would  most  likely  result  in 

A 

greatly  increased  computation  requirements.     The  second   is  how  R  varies 
with  H.     The  quantity   N[trace  (R)]    is   actually  2J(Z,U:H,R)   when  R  is  set 

A 

equal  to  the  identity  matrix.   As  seen  in  Appendix  D,  H  is  the  value  of 
H  such  that  the  matrix  gradient  of  J  with  respect  to  H  is  zero.  Then 

N  jjj  [trace  (R)]    =     -^  J2J(Z,U:H,R)/        1    =  0 

A  A 

or   in  some  respect  R  is   not  sensitive  to  small  errors   in  H. 
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4.3      Least  Squares    Identification 

The  Gaussian  distribution  of  v(n)   was    assumed   in  order  to  provide 
a  probalistic   framework,   however  the  same   equations   can  be   obtained  by- 
assuming   no  probability  model   for  v(n)    and    instead  relying  on  the 
nearly  universal  method   of  weighted   least  squares.      The   sum  of  weighted 
square   errors  between  z(n)    and  K  u(n)    is 

,N 


S    =  i    2    {zM    -  H  u(n)]'   c(n)  T      [z(n)    -  H  u(n)] 
n=l 


where    _T_     is   an  unknown  weighting  matrix.     Since  S   is  the  same  as 
J(Z,U:H,R)  when     T   =  R~    ,   from  Appendix  D  it  can  be  seen  that  the 
estimate    for  H  did   not  depend   on  R  anyway  and   the  value   of  H  which 
minimizes  S    for  given  sets   z  (l)  ,z(2)  , .  . .  ,z(N)    and  u(l)  ,u(2), .  . .  ,u(N) 
will  be  the  same  as   the  maximum  likelihood  estimate   for  H.      Under  such 
circumstances,   a  reasonable  estimate   for  E^c  (n)v(n)v' (n)v  is 

N 
n=l 

=  1   =1   2    c(n)[z(n)    -  Hu(n)][z(n)    -  H  u(n)] ' 
n=l 

which  is  the  same  as  the  maximum  likelihood  estimate  equation  (4-.3)« 


R  =E{c(n)v(n)v'(rO)  *  \   2,  c  (")v(n)0»  (n) 
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L.U     Error  Analysis 

The  error   analysis  which   follows   only  uses   the   covariance  of  v(n). 
For  both  the   approach  assuming  Gaussian  noise   and  the   least   squares 
approach,   the   covariance   of  v(n)  was    assumed  to  be  E(v(n)v'(n)v  =      r  r  R. 
Therefore  results   of  the   error   analysis   are   applicable  to  both   approaches 

When  the   equation  for   z(n)    is   substituted   into  equation   (4-. 2) 


N 


H   =  I   2   CMH  u(n)u'(n)  + 
n=l 


I   "  H   =      £    c(n)v(n)u'(n) 
n=l 


n  -\r  n  i  -i 

S    c(n)v(n)u'(n)       2T   c(n)u(n)u»  (n) 
n=l  ]ln=l  J 

'  N  1-1 

5,  c(n)u(n)u'(n)l  U.i 

n=l  J 


One  of  the  conditions  on  the  solution  was  that  v(n)  was  not  correlated 
with  u(n).  It  is  interesting  to  determine  the  effect  when  that  condi- 
tion is    not  met.      Let  W  be   a  weighted   average   of  u(n)u'(n),    n  =  1,2,...,N 


m        -i-i  r  n  "\ 

W  =      2    c(")  X   c(n)u(n)u'(n) 

.n=l  J         U^l  J 


U.5) 


If  v(n)  is  correlated  with  u(n)  so  that  the  weighted  average  of 
v(n)u'(n),  n  =  1,2,. . .  ,  N,  is 


M 


2  cm 

-1     " 

n=l           J 

N 


n=. 


(n)v(n)u'(n)|  ^  z£v(n)u' (n)^ 


Then 


H  -  H  =  M  W 


«-l 


(4.6) 
U.7) 


and  as  N  increases  H  -  H  will  not  go  to  zero  as  would  be  the  case  if 
v(n)  were  not  correlated  with  u(n).  Error  is  introduced  in  the  form 
of  a  bias  which  is  proportional  to  the  correlation  of  v(n)  with  u(n) 
and  inversely  proportional  to  the  "power"  of  u(n). 

From  equation  (£. 2)    it  may  be  noted  that  only  the  i   element  of 


z(n)  affects  the  estimate  of  the  i  '  element  of  H.   Similarly  from 

+  Vi  tin 

equation  (l.l)   only  the  i   element  of  v(n)  affects  the  i   row  of 

A 

the  error  H  -  H.   This  reflects  the  fact  that  K  was  assumed  to  be, 

in  general,  not  symmetric,  and  that  the  i   element  of  z(n)  is  affected 

only  by  the  i   row  of  H  and  the  i  element  of  v(n).   Let  h.  be  the 
ith  row  of  H. 


H  = 


t 
*2 


*K 


1L 

Then  the  estimate  of  the  i   row  of  H  is 


A  t 

h.  = 
— 1 


^  c(n)z. (n)u' (n) 
n=l      1 


2  c(n)u(n)u'(n)l 


-1 


U.8) 


and  the  error  is 


K  1  I 

h.  -  h.  = 
-i   -i 


N 


£  c(n)v  (n)u'(n) 
n=l     x 


N 


c(n)u(n)u'  (n) 


-1 


U.9) 


th        .th       A 

The  error  covariance  between  elements  of  the  i   row  and  j   row  of  H 


is 
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=  E 


(    2  c(n)u(n)u'(n)]        jg    c  (n)u(n)v.(n)l 

2   cCnVj^u'di)    J    J  2   c(n)u(n)u'(n)l      | 

N  "i-l   f    N       N  I 

Z   c(n)u(n)u'(n)  £     2  c  (n)u(n)u'  (m)c  (m)E$v.  (n)v.  (m)M 

n=l  J       Ln=l  m=l  t  l         J       >j 

"    N  -[-1 

x        2    c(n)u(n)u'(n) 
.n=l  J 

Since  Efv(n)v»  (m)j   =  0  for  n  /  a  and     Efv(n)v'  (n)7    =     -t-t  R 

-ij    =  E((-i   "  ~i)(^j    "  V'}   =  Rij   I  2   c(n)u(n)u'(n)J 

Let  c  be  the   average   of  the   c(n)    Ts. 

1        N 
n=l 

Then  in  terms   of  W,   the  weighted  average   of  u(n)u' (n)   defined  by  equation 

(A. 5)   the   error  covariance   is 


U.  io) 


R 


11  v"1 


-ij    N   c  - 


(A. 11) 


As  might  be  expected,  on  the  average,  the  error  covariance  of  the  esti- 
mate is  inversely  proportional  to  the  number  of  measurement  sets.   In 
terms  of  signal  to  noise  ratios,  the  quantity  c  R. .  is  a  measure  of 
the  pover  of  the  noise  v(n),  and  VJ   is  a  measure  of  the  power  of  the 
signal  u(n). 
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5.   IDENTIFICATION  OF  POWER  SYSTEMS 

The  problem  of  identifying  the  equivalent  admittance  matrix  and 
estimating  statistics  of  the  equivalent  bus  power  injection  vector 
for  the  external  system  model  has  been  placed  into  a  general  form, 
and  the  solution  to  the  general  form  has  been  given  subject  to  certain 
conditions.  The  two  models  for  identifying  the  equivalent  system  are 
The  Tie  Line  Power  Flow  Model 

PTL(n)  =  (YTL  +  YE^X)DB(n)  -  PEQX(n) 

The  Boundary  Bus    Impedance  Model 

D3(n)    =  ZBB[pEQl(n)  +   PB(nj^    +   ZB3   PEOJC(n) 

where   n  =  1,2 , . . . ,  N. 

In  both  models  variables  such  as  DB(n)    are  the  changes   which  occur 

between  times   t      and   t      .,  !B(n)   =DB(t    ,,)    -D3(t   ).      The  general 
n  n+1     —  —     n+1         —     n  E 

problem  is 

z(n)   =  H  u(n)  +  v(n)  n  =  1,2,...,N 

E{v(n)v.(n)}=     ^y  R 

and  the  solution 

N  ir-N  1-1 

21  c(n)z(n)u'(n)       £  c (n)u(n)u» (n) 
[n=l  iln=l  J 

N  • 

R  =    |  2  c(n)[z(n)   -  H  u(n)][z(n)   -  H  u(n)] 


A 

H  = 


n=. 


was  based   on  the   assumptions 
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1.  Elements   of  H   are   independent   of  each  other. 

2.  E^v(n)^   =  0 

3.  E^v(n)  ▼«fon)j   =  0     for     n  /  m 

L.  E{v(n)v'(n)]=     ^y  R  n=l,?,...,N 

5.  E^v(n)  u'(n)j   =  0 

Now  the  effects  of  these  assumptions  on  each  model  vill  be  discussed. 
For  both  models  if  symmetry  of  H  were  taken  into  account  at  the 
beginning  it  would  be  necessary  to  perform  matrix  inversions  of  dimen- 
sion -  K(K+l)  where  K  is  the  dimension  of  H.   In  developing  the  error 

A 

covariance  equation  for  K  it  was  noted  that  with  K  not  symmetric,  each 
row  of  H  is  actually  identified  independently  of  the  other  rows . of  H. 
With  symmetry  this  is  no  longer  true  and  the  result  is  that  all  inde- 
pendent elements  of  H  (either  upper  or  lower  diagonal  must  be  identified 
together.  The  size  of  the  matrix  inversion  necessary  can  be  significant 
as  the  dimension  of  H  increases.   For  example,  if  the  dimension  of  H 
is  10,  the  dimension  of  the  matrix  inversion  is  10  when  H  is  not 
assumed  to  be  symmetric  because  8S  was  shown  earlier,  each  row  of  H  is 
identified  independently  of  the  other  ro^s  and  only  one  matrix  is 
inverted.   However,  an  inversion  of  dimension  55  is  required  when 
symmetry  is  accounted  for  since  the  55  independent  elements  must  be 
identified  together.   One  way  to  handle  symmetry  at  the  end  is  to 
treat  the  ij   element  and  the  ji   element  of  H  as  two  estimates  of 
the  ij  '  element  of  H  (for  i  /  j).   Then  using  the  corresponding  error 

variances,  the  weighted  estimates  can  be  combined.   Let  a  be  the  ij 

A  th  A 

element  of  ][  and  r  its  error  variance,  and  b  the  ji   element  of  H 
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and  q  its  error  variance,  then  the  combined  estimate  c  is 

a  b 
-  +  - 
r     q 


c  = 


^  ♦  I 
r     q 


Symmetry  of  ZBB  and  YEQX  is  merely  the  result  of  reciprocity  since 
they  represent  passive  networks  (the  sources,  FEQI,  F3 1  and  PEQX  have 
been  moved  outside).   However,  in  the  case  of  ZBB  what  is  ultimately 
important  is  the  product  ZBB  (FEQI  +  P3_).   By  the  principle  of  recip- 
rocity an  input  at  terminal  1  causes  an  output  at  terminal  2  equal  to  the 
output  at  terminal  1  caused  by  the  same  input  at  terminal  2.  But  because 
of  different  demands  throughout  the  network,  the  input  at  terminal  2  may 
never  be  as  large  as  the  input  at  terminal  1. 

Had  the  actual  power  injection  and  voltage  angles  been  used,  v(t  ) 
(i.e.  PEQX (t  )  or  ZBB  PEQX (t  ))   would  most  likely  have  a  non  zero  mean 
value.  However,  because  changes  in  values  were  used,  it  may  be  expected 
that  v(n)  will  have  a  zero  mean  value  provided  all  the  measurements  are 
not  taken  when  the  entire  system  is  moving  in  the  same  direction.   For 
the  case  when  the  entire  system  is  moving  in  the  same  direction,  it  may 
turn  out  that  E<v(n)v' (n-k)V  =  0  for  k  greater  than  1  or  2.  That  is, 
v(n)  is  only  correlated  with  only  the  last  one  or  two  measurement  sets. 
This  could  be  accounted  for  by  a  modification  to  the  approach,  but  it 
is  not  anticipated  that  time  correlation  of  v(n)  will  be  a  major  source 
of  error.   Also,  when  the  entire  system  is  moving  in  one  general  direc- 
tion, depending  on  the  sampling  interval,  there  will  still  probably  be 
some  independent  movement  of  certain  buses  in  the  opposite  direction 
due  to  base  loaded  generators,  load  fluctuations,  etc. 
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The  form  of  R(n)  =  E^v(n)v»  (nft  =  — 7~y  R  was  based  on  the 

assumption  that  the  covariance  of  the  noise  is  dependent  only  on  the 

time  interval  between  measurement  sets  z(t  ),  u(t  )  and  z(t  , )»  u(t  .). 

—     n  '  -     n  -     n+1"  -     n+1' 

It   is  then  reasonable  to   assume  that   the  covariance   of  v(n)    =  v(t      .)-v(t   ) 

-  -     n+1     -     n 

is   of  the   form     R(n)   =   (t     ,    -  t   )R.     There  may  be  error  in  the  measure- 

-       n+1    n  —  J 

ments  z(n)  and  u(n)  but  it  was  shown  how  such  error  could  be  lumped  with 
v(n).   If  v  (n)  and  v  (n)  are  the  errors  of  z(n)  and  u(n)  respectively, 

Z  U 

then  z(n)    =  H  u(n)   +  \f(n) 

where     v(n)    =  -v   (n)  +  H  v   (n)  +  v(n) 
-        .       -z  u  - 

Since  v(n)    is   independent   of  v    (n)    and  v    (n),  E< v(n)v' (n)V    =  2  +      "  7    V  5 
where     2  =  E  Hi  I   (n)    -  Y.   M]   [M  v   (n)    -  v   (n)]    Y     and  should  be 
relatively  constant  since   it  depends   on  measurement  error.     The  measure- 
ment error  should  be  in  the  vicinity  of  one  to  ten  percent  while  v(n) 
can  be  expected  to  be   in  the  vicinity  of  IGOJf  of  H  u(n)    so  the  effect  of 
2  on  the   estimate  of  R  should  be  very  minor. 

The   last   assumption,   that  v(n)    is   uncorrelated  with  u(n)   can  be   a 
source   of  considerable   error  when  this   condition  is    not  met.      The  error 

A 

between  K  and  H  was   shown  to  be 


r  n 

H   -  H   =    MT    c(n)v(n)u'(n) 


'  N  -i-l 

2    c(n)u(n)u»(n) 
n=l  J 


For  the  boundary  bus    impedance  model 

DB(n)    =  ZB3[FEQl(n)  +   P3(n)]    +   Z3B   FEQX(n) 

v(n)u'(n)    =  ZBB   PEQX (n)[lSQl(n)  +   P3(n)~) 

=  ZB3  AQX[PX(n)PIf(n)AQI'  +   PX(n)F3'(n)] 

while  for  the  tie  line  power  flow  model 
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FTL(n)    =    (YTL  +  YEOX)PB(n)    -  PEQX(n) 

v(n)u»(n)    =   -PEQX(n)D3» (n) 

=  -AQX  RL(n)[AQI  Pl(n)   +  FB(n)  +  AOX  FX(n)]'    ZBB 
=  -AQX[PX(n)PI'(n)AQP  +  PX(n)fB' (n)]  ZBB 

-  AOX  PX(n)PX' (n)AQX'ZPB 

where   AOJX   end   AOI  were   defined    in  equations    (3-7)   and    (3. 13).      Both 
models   have  the   term     AQX[PX(n)PT' (n)AOJ  +   PX  (n)FB  '  (n)"J      in  common. 
However,   the  tie   line   power  flow  model  has   the    additional  term 
AOX  PX(n)pX'(n)AOX'    =  PEQX(n)PEQX' (n).      Ordinarily   it  would  be   expected 
that   a  charge   in  power   injection  at   one   particular  bus  will  be   indepen- 
dent  of  the   changes    in  power   injections    at   all  other  buses   so  that 
s[pX(n)PI»(n)^   =  0     and     E[pX(n)F3 '  (n)}  =  0.     The  exceptions  might  occur 
when  the   entire  system   is  moving   in  the   same  direction  as   discussed 
earlier.      However   in  the  case  of  the  tie   line  power  flow  model, 
E|FE0X(n)PZ0X'  (nH   =  R     is   not   zero  unless   the   external  system  bus 
injections   do   not   change.      This   causes    a  bias    in  the   estimate   of  H  which 
depends   on  the  magnitude   of  FEQX(n)    in  relation  to  FEQl(iQ  +   FB(n),  but 
in  general   it  can  be   expected  that   the  two  will  be   of  comparable  magni- 
tude so  that  the  bias    introduced   is    in  the   vicinity  of  100/6.      Using  W 
and  M  as   defined  by  equations    (4-. 5)    and    (4. 6) 

r  n        i-i  r  n 

M  =  |2    c(n)l  2  c(n)v(n)u'(n) 


f£  l"1 

E    c(n) 
n=l  J 


N 

2    c(n)FE0X(n)PE0X'(n) 
n=l 


Z?B 


«s-R   ZBB 
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the  bias  in  the  estimate  for  the  tie  line  power  flow  model  is 

H  -  H  =  M  W   &-R   ZBB  W"1 

A         A 

so  that  the  estimation  equations  for  H  and  R  are  really  coupled  for 
this  model. 

The  boundary  bus  impedance  model  and  the  tie  line  power  flow  model 
are  two  specific  approaches  to  the  problem  of  identifying  the  equivalent 
external  system  model.   In  the  former  the  inputs  are  Own  System  bus 
power  injections  and  the  outputs  are  boundary  bus  voltage  angles.   In 
the  latter  the  inputs  are  boundary  bus  voltage  angles  while  outputs  are 
tie  line  power  flows  (which  depends  on  voltage  angles).   Because  it  is 
power  which  is  bought  and  sold,  while  voltage  angles  are  determined  by 
the  line  admittances  and  the  distribution  of  bus  power  injections,  it 
can  be  seen  that  bus  power  injections  are  the  legitimate  independent 
inputs  and  not  voltage  angles.  This  is  the  cause  of  the  bias  in  the 
estimate  for  the  tie  line  power  flow  model.  Both  models  need  D3(n)  but 
the  tie  line  power  flow  model  uses  the  tie  line  power  flow  vector, 
PTL(n),  which  in  general  will  be  more  accurate  a  measurement  than  the 
estimate  for  the  equivalent  injections  from  IS,  FEQl(n) .   It  is  possible 
to  measure  PTL(n)  directly  while  the  actual  value  of  PEOl(n)  is  dependent 
on  voltage  magnitudes  and  transmission  line  resistance  as  shown  by  the 

A        A 

nonlinear  analysis  in  Appendix  C.   However  the  solutions  for  H  and  R 
are  more  readily  found  for  the  boundary  bus  impedance  model. 

Obtaining  the  estimate  of  YE OX  from  the  estimate  of  ZBB  is  straight- 
forward.  However,  if  the  estimate  of  YE  OX  is  to  ultimately  be  used  in 
a  linear  load  flow  model  for  OS,  then  as  shown  in  Appendix  B,  the 
estimate  for  ZBB  can  be  used  directly.   If  ZOS  is  the  bus  impedance 
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rcatrix  for  OS, 

DI 

DB 


zos 


F3  +   FEOX 


where 


ZOS   = 


YII"1       0 


P.OV 


ZB3  Iaoj      f) 
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6 .   S  IMULAT 10  N  RES  ULTS 


In  order  to  test  the  method  proposed,  a  computer  simulation 

program  was  developed  to  generate  data,  identify  the  equivalent  model, 

and  check  how  well  the  identified  model  could  predict  line  flow 

changes  in  OS.   Figure  6.1  gives  a  broad  overview  of  how  the  simulation 

was  conducted.   First,  using  transmission  line  data,  matrices  to  be 

used  in  the  simulation  were  coirputed.   These  include  ZW§  (impedance 

matrix  for  Whole  System),  AQI.  YIl"1,  (YEQI  +  TBOS),  (YEQX  +  YTL),  AQX. 

and  YTL.  For  convenience,  the  confidence  coefficients,  c(l),  c (?),..., 

c(N),  were  all  unity.   Matrices  CHI  and  SGIf-V  accumulated  a  running  sum 

of  z(n)u'(n)  and  u(n)u'(n)  respectively.   In  order  to  simulate  the 

variations  in  bus  power  injections  throughout  the  network,  P(n)  (vector 

of  bus  power  injection  changes  for  WS)  was  generated  by  random  numbers 

having  an  RMS  (root  mean  square)  value  which  was  a  fraction  of  the 

nominal  operating  value  of  bus  power  injections  P   .   One  fraction, 
1  ^  r  tt  -nom 

FCTOS,  was  used  for  all  buses  in  OS  and  another,  FGTXS ,  was  used  for 

all  buses  in  XS.   (i.e.  E  «  EX?(n)l  =  fpCTXS  x  FX    ~|2)   This  main- 

C      i       J  L  nom . J 

l 

tained    a  relative  scale   for  changes  while  at  the  same  time   it  simulated 

independent   changes    in  power   injections.  As   shown  in  Section  3>   when 

the   linear  model  is   used,   the  bus  voltage  angle  vector   for   OS   can  be 
found   using   either 


Dl(n) 
D3(n) 
DX(n) 


ffiJS 


Pl(n) 

F3(n) 
PX(n) 


(6.1) 


or  equivalent ly 
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Bl(n) 


r^ 


3(n) 


'fl(n) 


IB'   Y30S 


FB(n)  -  PTL(n) 


(6.2) 


Ordinarily  equation  (6.2)  or  its  nonlinear  analog  would  be  used  to 
estimate  the  state  of  OS  using  real  measurements  of  Pl(n),  P3(n),  and 
PTL(n)   (and  any  other  redundant  measurements).   However,  for  purposes 
of  simulation  equation  (6.1)  would  have  to  be  solved  first  to  find 
DX(n)  in  order  to  generate  PTL(n).  Since  either  method  results  in 
the  same  values  of  D3(n)  (in  the  linear  model),  equation  (6.1)  was 
used  to  simulate  the  output  of  a  state  estimation  for  ]2B(n).  The 
simulated  state  estimation  for  OS  results  in  errorless  data.   However 
as  discussed  earlier  the  only  effect  of  measurement  noise  and  state 
estimation  errors  is  to  modify  the  disturbance  vector  v(n). 

Everything  described  up  to  now  was  for  the  purpose  of  generating 
data  for  the  identification  routine.   Using  generated  data,  Pl(n) , 
PB(n),  and  D3(n),  the  input  and  output  measurement  vectors  were  formed, 
u(n)  =  A^I  Pl(n)  +  P8(n)  and  z(n)  =  DB(n).  These  values  were  also 
stored  for  later  use  in  finding  R.   Running  sums  of  z(n)u'(n)  and 

u(n)ur(n)  were  maintained  by  CHI  and  SGI"/  so  that  at  any  time  the 

A  -1         A 

current  estimate  H  could  be  found  by  H  =  CHI  SGIW   .   Using  H  and  the 

stored  measurement  sets,  an  estimate  K  was  found.   Parallel  to  all  of 

this,  the  actual  disturbance  vector,  v(n),  found  by  v(n)  =  ZB3  FEQX(n) 

=  ZBX  PX(n),  was  used  to  maintain  a  running  sum  of  v(n)v'(n)  in  R. 

The  estimated  values  H  and  R  could  then  be  compared  to  the  actual 

values  H  =  Z33  and  R.   In  order  to  compare  how  well  the  identified 

model  could  predict  changes  in  line  power  flows,  the  estimate  for  ZOS , 
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the  impedance  matrix  relating  OS  bus  power  changes  to  OS  bus  voltage 
angle  changes  vas  found  by 


zos  = 


YII 


0 


-1 


0 


/.or 


ZB3  [AQI   I] 


(6.3) 


Again  random  variations  of  bus  power  injections  were  generated  using  a 
nominal  base  value  and  estimated  bus  voltage  angle  changes  found  by 


lk(n) 


D3(n) 


=  ZPS 


PI(n) 

F3(n) 


(6.0 


were  used  to  estimate  line  flow  changes.   The  actual  bus  voltage  angle 
changes  using  the  same  bus  power  changes  were  found  from  D(n)  =  Z1-/S  P(n) 
and  were  used  to  find  actual  line  flow  changes  which  could  be  compared 
to  the  values  found  using  the  identified  model. 

Throughout  the  program,  in  order  to  present  data  in  a  condensed 
form,  for  many  variables  (z(n),  u(n),  v(n),  actual  line  flows,  errors 
in  predicted  line  flows,  and  generated  bus  power  changes)  only  the  RNS 
value,  the  mean,  the  minimum,  and  the  maximum  values  were  maintained. 
It  should  also  be  pointed  out  that  the  basic  approach  is  simple.   Only 
a  small  portion  of  the  program  was  involved  in  the  actual  identification 

A         A 

of  H  and  R.  The  rest  is  support  for  generating  data  and  evaluating 
results. 

Using  equation  (£. 11) i  the  estimated  error  covariance  matrix,  P. . , 
of  row  i  of  H  is  for  c(n)  =1.0 


K  1  *      -1 

P   =  -  R   W 
-ii    N  ii  - 
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where 


W 


-1 


\l  t    u(n)u'(n)l 
L     n=l  J 


-1 


=   Nx  SGIMT"1  =  Nx  S2GMA 


Hence  the   estimated  standard  deviation  of  element   ij   of  H   is 

(R..   x  SIGMA..)    .      Since  R.  .   should  be  somewhat  proportional  to  the 
11  jj  n  v      v 

square   of  PCTXS  while  SIGMA.,   should  be   inversely  proportional  to  the 

J  J 

square  of  PCTOS,  in  the  case  of  the  linear  model  with  perfect  measure- 
ments of  z(n)  and  u(n),  identification  error  is  dependent  on  the  ratio 
of  PCTXS  to  PCTOS  (percentage  variation  of  XS  bus  powers  to  percentage 
variation  of  OS  bus   powers). 

One  system  presented   in  a  paper    C  2  J  was   used  here  to    test  the 
method.      The  system  shown  in  Figure    (6.2)   consists   of  18  buses.      For 
purposes   of  simulation  Whole  System  consists   of  these   18  buses   of  which 
buses   l,2,3j£,    and    18  belong  to   IS ,    5  through  7  belong  to  BS ,    and  8 
through   17  belong  to  XS.      Bus   18   is  the   reference  bus.      Other  parti- 
tioning  could  be  used  but   this   one   does  result   in  matrices    (YEQX  +  YTL) 
and    (TEQ I  +  YJ30S)     which   are   comparable   in  size.      This   avoids    a  case   in 
which  coupling  of  BS   with    IS    is  very  much  stronger  than  BS   with  XS. 
For  this   system    (YEQI  +  YBOS)    and    (YEQX  +  YTL)    are 


YEOI  +  YBOS    = 


13.63       0.0  -A. 07 

0.0         9.^3         -9. A3 
-U.  07     -9./+3         85.9 


YEQX  +  YTL     = 


11.5       -6.M        -5.36 
-6.H     H.3  -8.12 

-5.36     -8.12         13.5 
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Both  the  effects  of  increasing  disturbances  and  increasing 
measurements  are  simulated  here.   In  Figure  6.3  and  Table  6.1  are  the 
results  of  increasing  disturbances  (PX(n)).  The  RMS  of  bus  power 
changes  for  OS  in  each  case  was  10?  of  nominal  while  RKS  values  of 
PX(n)  were  Of.  9    1.25?,  2.5*,  5.0g,  1056,  and  ?C?   of  PX   .   Another  way 
to  say  this  is  that  "signal  to  noise  ratios"  were  co  :1,  3:1,  4*1»  2:1, 
1:1,  and  1:2  respectively.   The  same  sequence  of  random  numbers  were 
generated  each  time.   Included  in  Table  6.1  are  actual  values  of  ZBB, 
actual  errors  in  identifying  it,  and  the  estimated  standard  deviations 
using  R  in  equation  (4..  11).   Also  included  are  the  values  of 
(YE  )X  +  YTL)  estimated  from  ZBB,  and  the  errors  of  predicting  changes  in 
line  flows  using  ZBB.   In  a  similar  format,  the  results  of  increasing 
the  number  of  measurements  are  summarized  in  Figure  6.4  and  Table  6.2. 
Changes  in  both  OS  and  XS  bus  power  injections  had  RKS  values  of  10? 
of  the  nominal  operating  point  while  H  was  identified  using  16,  32,  64, 
12s?,  256,  and  512  measurements. 

The  results  show  that  not  only  is  identification  error  commensurate 
with  the  noise  and  number  of  measurements,  but  also  estimated  standard 
deviation  gives  an  accurate  measure  of  the  identification  error. 
Because  the  same  random  numbers  were  generated  for  each  value  of  PCTXS 
in  the  first  case,  the  actual  errors  shown  in  Figure  6.3  are  nearly 
exactly  linear  since  as  was  shown  earlier,  for  equal  numbers  of 
measurements  identification  error  depends  on  the  ratio  of  FCTXS  to 
PCTOS.   For  the  case  of  increasing  measurements  in  Figure  6.4,  ZBB-.. 
does  not  have  a  steadily  decreasing  error  as  with  ZBB   ,  but  it  is  still 
in  the  range  of  the  estimated  standard  deviation. 
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To  contrast  identification  models,  in  Table  6.3  the  results  of 
identifying  YE  OX  +  YTL  by  the  tie  line  power  flow  model  along  with  the 
estimated  standard  deviations  are  compared  with  YEQX  +  YTL  as  found 
from  ZB3  which  was  identified  by  the  boundary  bus  impedance  model 
using  the  same  data.   Besides  identifying  the  wrong  values,  the  tie 
line  power  flow  model  gives  inaccurate  estimates  of  error. 

In  Figure  6.5  is  another  network  presented  in  a  paper  by  Stagg  C9  J. 
OS  consists  of  the  higher  voltage  network  and  XS  is  the  lower  voltage 
network.  The  "tie  lines"  are  really  transformers.  Using  the  boundary 
bus  impedance  model  an  equivalent  network  is  identified  to  replace  the 
21  buses  of  XS,  so  that  ^S  can  be  analyzed  with  9  buses  instead  of  the 
original  30.  The  results  in  Table  6.£  again  show  close  correspondence 
between  actual  error  and  estimated  standard  deviation. 

The  simulation  carried  out  was  certainly  not  fully  realistic. 
However  the  realistic  aspects  include 

1.  Values  of  bus  power  changes  were  found  first,  then  a  state 
estimation  for  OS  was  simulated  and  appropriate  values  were 
passed  to  the  identification  portion  of  the  program. 

2.  The  effects  of  increasing  measurements  were  studied  under 
the  assumption  that  XS  would  vary  as  much  as  OS. 

The  simulation  was  unrealistic  in  that 

1.  The  nonlinear  model  was  not  used  to  generate  data  with  the 
result  that  a  linear  identification  model  was  used  to  find 
parameters  to  a  linear  system.  To  generate  data  from  a 
nonlinear  model  properly  would  require  a  complete  load  flow 
for  each  measurement  set.   (it  would  be  simple  to  vary  voltage 
magnitude  and  voltage  angle  and  then  find  bus  power,  but  bus 
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power  must  be  the  specified  variable  as  it  is  the  independent 

input,  not  voltage  angle.) 
2.  Bias  power  injections  were  independently  generated  whereas 
realistically  they  will  probably  be  correlated  in  space  and 
time.   There  were  no  conditions  on  the  solution  for  correla- 
tion in  space.   In  fact  FEQX(n)  itself  has  a  covariance  matrix 
E  ^FEQX (n)PEQX ' (n) r    which  is  in  general  not  diagonal.   However, 
the  solution  did  assume  no  correlation  in  time.  P.s   discussed 
earlier  in  Section  5,  this  effect  can  be  accounted  for  by 
modifying  the  approach. 
Based  on  the  results  of  these  and  other  simulations  certain  con- 
clusions can  be  made  about  the  boundary  bus  impedance  model.   First, 
with  enough  measurements  it  is  possible  to  identify  the  equivalent 
model  under  conditions  of  the  simulation  (independent  movement  of  the 
bus  injections).   The  number  of  measurements  may  seem  rather  large  but 
of  course  the  noise  is  also  rather  large.  The  results  show  that  as 
expected  either  a  decrease  in  the  disturbances  by  a  half  or  using  four 
times  as  many  measurements  approximately  halves  the  error.   Second  and 
possibly  even  more  important,  the  error  predicted  by  using  the  estimated 
R  in  equation  (£.11)  agrees  very  well  with  the  actual  answer.  This 
means  that  an  accurate  estimate  can  be  predicted  from  only  the  rceasure- 
ments  used  to  identify  H  so  that  the  limits  to  which  H  can  be  trusted  are 
also  known.  This  is  of  course  not  the  case  with  the  tie  line  power  flow 
model. 

A  very  important  question  which  may  arise  is  how  rr.uch  data  would 
be  required  to  identify  the  system.  The  simulation  was  carried  out 
under  idealized  circumstances  and  as  was  shown  earlier,  for  the 
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simulation  the  identification  error  for  equal  numbers  of  measurements 
depends  only  on  the  ratio  of  XS  bus  power  changes  to  03  bus  power 
changes.   In  actuality  however  there  are  upper  and  lower  limits.   If 
data  were  used  from  samples  taken  too  close  together  in  time  so  that 
changes  were  small,  the  effect  of  measurement  noise  and  state  estima- 
tion error  could  become  very  significant.   On  the  other  hand,  the 
system  can  vary  over  only  so  much  of  a  range,  so  that  if  the  time 
between  measurements  were  too  large,  then  because  of  cyclical  patterns 
of  power  requirements,  there  might  not  be  much  new  information  after 
the  first  few  measurements.  To  determine  the  amount  of  data  required 
and  the  best  interval  between  samples  would  require  a  study  using 
actual  operating  data. 
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y\ 

YEQX  +  YTL 

Actual  Value 

11.50 

-6.14 

-5.36 

Error  BBIM 

.CO 

-1.37 

1.30 

Error  TLPFM 

.30 

-14.44 

10.66 

E.S.D.  TLPFM 

.19 

.40 

.39 

Actual  Value 

-6.14 

14.26 

-8.12 

Error  EBIM 

.07 

-1.87 

2.97 

Error  TLPFM 

4.52 

-19.37 

14.18 

E.S.D.  TLPFM 

.24 

.52 

.51 

Actual  Value 

-5.36 

-8.12 

13.47 

Error  BBIM 

.52 

-  .52 

2.79 

Error  TLPFM 

5.41 

-23.24 

16.07 

E.S.D.   TLPFM 

.60 

1.29 

1.26 

Actual  Value   —  Actual  Value   of  YEQX  +  YTL 

Error  3B TM       --  Error   in   Identifying  YEQX  +  YTL  using  the  Boundary 

Bus    Bapedance  Model  YEQX  +  YTL  =  ZBB   -    (YEQI  +  Y30S) 

Error  TLPFM     --  Error  in   Identifying  YEQX  +  YTL  using  the  Tie    Line 
Power  Flow  Model 

E.S.D.  TLPFM  --  Estimated  Standard  Deviation  using  R  from  the  Tie   Line 
Power  Flow  Model 

OS   and  XS  varied  equal  amounts    (10^   of  nominal)    and   the  same 
256  measurements  used   for  both  models. 


Table  6.3     Comparison  of   Identification  Models 


Z33 

for  30 

72- 

3  us 

System 

Actual  Value 
Actual  Error 
E.S.D. 

.0644 

-.0016 

.0019 

.0489 
.0004 
.0020 

.0491 
.0179 
.0217 

Actual  Value 
Actual  Error 
E.S.D. 

.0#?9 

-.0019 
.0021 

.0665 
.0007 
.0023 

.0663 
.0184 

.0244 

Actual  Value 
Actual  Error 
S  S  D 

.0491 

-.0021 

.CC25 

.C663 
.0010 
.0026 

.1122 
.0222 
.C2?2 

E.S.D.    —  Estimated  Standard  Deviation  using  R. 
Z3B    Identified  "by  the  Ecundary  Eus    Impedance   --'odel 
lCr-  variation  of  bus   pcvers    in  both  XS   and   OS. 
2$6  measurement  sets  used. 


Table  6.4     Identification  of  a  Lover  Voltasre    Tetvork 
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7.   USE  IK  POWER  SYSTEMS— SOKE  PRACTICAL  ASPECTS 

7.1  Use  of  Bios  Power  Inputs 

From  equation  (£.10)  which  relates  the  error  covariarce  of  the 
estimate,  it  is  obvious  that  the  input  directly  affects  the  error. 
For  the  boundary  bus  impedance  model  u(n)  =  PEQl(n)  +  PB(n).   As 
discussed  previously  and  shown  in  Appendix  C,  FE3l(n)  really  depends 
on  the  operating  point  in  a  nonlinear  manner  so  that  it  is  possible  for 
FE3l(n)  to  introduce  significant  error.   The  error  from  P3(n)  on  the 
other  hand  is  only  due  to  measurement  noise.   For  this  reason,  it 
should  turn  out  that  the  identification  method  presented  will  be  most 
successful  when  the  boundary  buses  have  relatively  large  power  injec- 
tions (generators  or  loads)  which  supply  much  of  the  known  variation. 
Or  in  other  words,  it  would  be  nice  to  have  good,  strong  input  signals. 
In  actuality  IS  may  be  of  relatively  large  size  so  that  contributions 
of  Pl(n)  to  PEQl(n)  =  API  PT(n)  from  far  across  the  network  may  be  of 
questionable  value  and  accuracy  due  to  the  linear  approximation  of  API. 
Since  intuitively,  buses  in  IS  closest  to  BS  will  have  the  most  affect 
on  PEOl(n)  and  involve  the  least  error,  it  may  be  practical  to  use  only 
part  of  the  contribution  of  Pl(n)  to  PEQl(n)  and  lump  the  rest  with 
PEQX(n). 


7.2  Vodel  Verification 

Although  contrary  arguments  could  be  made,  the  boundary  bus 
impedance  model  is  probably  the  better  of  the  tv.ro  models  since  it  does 
not  involve  estimating  a  significant  bias.   However,  the  tie  line  power 
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flow  model  is  still  valid  for  verification  purposes.   For  the  boundary 
bus  impedance  model  from  equation  (A. 3) 

i      N 
R=     "N2    c(n)0(n)0'(n) 
n=l 

=     tt   g   c(n)ZBB   FEQX(n)FEQX'(n)ZBB  (7.1) 

n=l 

A  ^-1 

If  R   is   premultiplied   and  post  multiplied  by  ZB3        then  an  estimate  of 
the   covariance   of  PS OX (n)    is 

HpEQX   =  ZB3     2  ZB3  (7.2) 

fbw  if  measurements    of  PTL(n)  were  taken  at  the   seme  times   as  measure- 
ments   of  D3(n)   and   FEOl(n)  +   F3(n),    and   if  the   estimate   of  YEQX   is 

YEQX   =  Z33   x   -  YEQI  -  YB3 

/\  a 

then  equation    (7.3)   evaluated   for  YEQX  should  be  close  to  Eptrny   °f  ■ 
equation    (7.2). 

1      N 

i   2  c(n)FEQX(n)PEOX'(n) 


N 

n=l 


N 


=     fr   T<  {[PTLCnl    -    (YTL+  YEQX)DB(n)] 

n=l   L 

x  [?TL(n)    -    (YTL  +  YEQX)D3(n)V    c(n)V  (7.3) 

If  the  two  were    not   close,   then  the  tie   line  power   flow  model  would  be 
indicating   a   lack  of  validity  for  the   estimate   of  ZBB   from  the  boundary 
bus    impedance  model.      The   tie   line  power   flow  model  could   also  be  used 
"on-line"   once  YEQX  had  been  identified  to  provide   a  continuing  check  of 
the  validity  of  the    identified  equivalent  system.     Should   the  measured 
value   of  FT L ( n)   vary   from    (YTL  +  YEQX)SB  (n)   an  amount   not  commensurate 
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1       ft 

with        /    c  Rp_„„,   that  would   indicate   a  change   in  XS  had   taken  place, 

The   change   of  course   could   either  be    in  transmission  line   status 
affecting  YEQX  or  a  significant  change   in  PX(n)   affecting  FEQX(n). 


7.  3 Estimating  Equivalent  Power    Infections 

Once  YEQX  has  been  identified,   the    actual  values   of  FEQX  (t   )   could 
be  estimated.     With  enough  measurements    it  might  be  possible  to  develop 
an  approximate  daily  pattern  for  FEQX  (Q   so  that  they  could  be  used  as 
pseudo  measurements   [  8  3  for  the  purpose   of  predicting   a   future  power 
flow  situation.      If  such  were  the   case,   the   equivalent  system  could  be 
extended   in  use   from  predicting  changes    in  OS   during  which  the   change 
FEQX(n)  would  be  small  to  use   in  predicting   actual  line   power   flows   for 
given  power  injections   in  OS. 


7. L      Inputs  Which  Decrease    Identification  Error 

It  has  been  assumed  that  the   identification  would  be  primarily 
passive   in  nature    (i.e.   by   "listening"   to  the   system)   so  that   any  power 
changes   or  voltage   angle   changes   used  to   identify  the   equivalent  system 
would  be  due  to  changes   in  the  system's   load  demands.      It  was   also 
argued  that     ~~7~\  E  was    a  reasonable   approximation  to  the    noise 
covariance  where  c(n)   is    inversely  proportional  to  the  time,  t     ,    -  t^, 
between  measurements.      However,    any  large,    relatively  fast   changes    in 
OS   would   provide   an   input,   output  measurement   set  with  relatively  little 
noise    (FEOX(n)).      These  changes  could   be   planned   or   accidental.      If 
planned   changes  were  used,    it  would  be  desirable  to  obtain  a  maximum  of 
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invorrr.ation  with   a  minimum  of  change.      For  this   case    it  might  be 

possible  to  formulate  the  problem  in  terms   of  minimizing   a  cost 

function  subject   to   a  constraint  such   as   finding  the   inputs 

u(n)      (n=l,2,...  ,N)   which  minimizes  the  cost 

c   =  fx  (error  in  H)  4   f  [u(n)u»  (n) ,    (n  =  1,2,...,N)} 

subject  to 

z(n)    =  H  u(n)  +  v(n)  (n  =  1,2,..  .,N) 

where    f.    and   f     are  scalar   functions. 


im 


7.5   Use  of  the  Equivalent  System  When  the  External  System  is  Known 

Although  the  linear  load  flow  model  was  used  to  solve  the  proble 
of  finding  equivalent  systems  with  no  prior  knowledge  of  the  external 
systems,  there  is  still  an  advantage  to  using  the  model  when  enough  is 
known  about  the  external  system  so  that,  for  instance  YEQX  and  PEQX 
could  be  found  analytically  for  a  particular  system.  The  linear  model 
clears  the  smoke  so  to  speak,  so  that  a  good  approximate  picture  can 
be  easily  obtained.   For  example  AQX  (and  AQI)  gives  a  general  overview 
of  how  power  is  routed  by  the  network  while  a  comparison  of  YEQX  with 
YB3  shows  at  a  glance  the  relative  coupling  between  the  two  networks. 
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7.6   Identification  of  the  Nonlinear  Kodel 

Use  of  the  linear  model  may  be  sufficient  to  predict  line  power 
flows  and  voltage  angles.  However,  in  order  to  predict  the  effect  of 
XS  on  voltage  magnitudes  within  OS,  the  full  nonlinear  model  is  neces- 
sary, developed  in  Appendix  C  are  the  equations  relating  complex  bus 
voltages  to  complex  bus  power  and  current  injections.   The  computation 
required  and  complexity  of  identifying  the  nonlinear  equivalent  model 
would  likely  be  very  much  more  than  for  the  linear  model.   The  increased 
accuracy  in  predicting  real  power  flows  would  most  likely  be  minimal, 
but  the  linear  model  cannot  predict  voltage  changes  due  to  the  struc- 
ture of  XS. 
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8.   CONCLUSIONS 

3.1  .Areas  for  Further  Study 

The  tie  line  power  flow  model  produces  a  biased  estimate  of  the 
equivalent  model  parameters,  however,  if  an  unbiased  estimate  could  be 
obtained  so  that  accuracy  was  comparable  to  the  boundary  bus  impedance 
model,  there  would  probably  be  an  overwhelming  advantage  to  using  the 
tie  line  power  flow  model.   One  of  the  most  important  reasons  is  that 
it  uses  variables  from  a  local  region  only.   FTL(t)  is  directly 
measurable  and  in  order  to  find  DB(t)  state  estimation  would  be  required 
only  in  a  small  area  around  Boundary  System.   (Power  coming  in  from  the 
rest  of  Own  System  could  be  treated  the  same  as  FTL(t)  for  purposes  of 
state  estimation.)   Also  the  linear  tie  line  power  flow  model  appears 
to  be  less  of  an  approximation  to  its  nonlinear  analog  (see  Appendix  C) 
than  the  linear  boundary  bus  impedance  model  is  to  its  nonlinear  analog. 
Therefore,  it  is  suggested  that  a  study  to  find  a  method  for  eliminating 
bias  in  the  tie  line  power  flow  model  would  be  worthwhile. 

The  simulation  carried  out  used  a  linear  load  flow  model  for  the 
purpose  of  generating  data  used  to  identify  the  equivalent  model.  This 
simplified  the  simulation  and  it  also  served  to  eliminate  errors  which 
would  be  caused  by  variations  from  ideal  conditions  (low  line  resistance 
and  voltage  magnitude  nearly  constant).   A  more  realistic  simulation 
should  be  made  to  find  how  data  from  a  nonlinear  model  affects  accuracy 
of  the  identified  linear  model  parameters.   Also  the  net  power  flow  summed 
over  all  the  tie  lir.es  is  usually  held  to  a  scheduled  value.   The  effect 
of  this  on  the  identification  should  also  be  studied. 

As  mentioned  previously,  the  identification  accuracy  could  be 
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increased  by  varying  the  generation  pattern  for  the  purpose  of  probing 
the  external  system.   However  it  would  be  necessary  to  maximize  informa- 
tion and  minimize  the  variation  required.   It  is  felt  that  one  step  in 
the  proper  direction  might  be  to  move  one  generator  at  a  time  since  this 
would  eliminate  any  interference  among  the  inputs.   A  look  at  the  boun- 
dary bus  impedance  model  written  in  terms  of  the  general  model 
z(n)  =  H  u(n)  +  v(n)  shows  that  if  elements  of  u(n)  were  all  zero  except 
for  element  i,  then  the  i   column  of  H  would  be  — j— r  Vz(n)  -  v(n)j. 
A  study  into  this  problem  should  result  in  a  more  sophisticated  solution 
to  a  rather  complex  problem. 

An  obvious  extension  of  identification  of  the  linear  model  is  the 
identification  of  the  full  nonlinear  model.   It  is  felt  that  the  com- 
plexities involved  would  be  an  order  of  magnitude  above  that  required 
for  identification  of  the  linear  model.   The  accuracy  gained  in  predicting 
voltage  phase  angle  changes  would  probably  be  small,  but  voltage  magnitude 
changes  which  have  been  ignored,  would  be  available  from  a  nonlinear 
model.   If  a  search  type  of  solution  were  usei,  then  of  course  the 
solution  to  the  linear  model  identification  would  give  a  starting  point. 


g.2  Summary 

The  problem  of  predicting  how  an  electric  power  network  will  interact 
with  other  networks  has  been  solved  by  reducing  the  entire  outside  world 
down  to  a  relatively  small  equivalent  model.  The  problem  is  complicated 
by  the  facts  that  (l.)  r.o  knowledge  of  parameters  or  variables  of  the 
outside  world  can  be  assumed,  (2.)  identification  must  be  primarily 
passive,  and  (3.)  "signal  to  noise  ratios''  can  be  expected  to  be  about  1:1. 
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A  linear  model  for  the  overall  system  was  developed  and  its  use 

justified  with  respect  to  the  more  conventional  nonlinear  model  by  con- 
sidering accuracy  versus  computation  requirements,  insight,  and  com- 
plexity of  the  identification  approach  required.  Two  models  were  then 
proposed  with  each  having  certain  advantages  and  shortcomings.  The  tie 
line  power  flow  model  uses  data  which  is  relatively  accurate  and  locally 
available,  but  it  produces  a  biased  estimate.   On  the  other  hand  the 
boundary  bus  impedance  model  produces  an  unbiased  estimate  but  it 
requires  r'ata  from  all  over  the  system  and  also  requires  manipulation 
of  the  equations  for  own  system.   It  was  argued  that  changes  which  occur 
in  power  and  bus  voltage  angle  between  two  times  should  be  used  with  the 
hope  of  obtaining  a  zero  mean  disturbance  and  also  to  avoid  the  use  of 
matrices  whose  elements  have  magnitudes  very  large  compared  to  the 
determinant  of  the  matrix.  Since  both  models  were  in  the  general  form 
of  a  linear,  static,  input-output  relation  with  additive  disturbances, 
a  solution  to  the  general  problem  was  found  under  assumed  conditions. 
Then  the  solution  was  applied  to  the  two  models  and  it  was  shown  how 
the  tie  line  power  flow  model  produces  a  biased  estimate.  The  simulation 
made  to  test  the  method  was  discussed  shewing  the  excellent  correspon- 
dence between  theory  and  simulation  results.   Finally,  practical  aspects 
of  the  problem  were  briefly  discussed  including  model  verification, 
estimation  of  the  disturbances,  system  probing,  and  use  of  the  linear 
model  reduction  method  for  purposes  other  than  identification. 

This  method  is  proposed  as  a  means  of  modeling  the  outside  world 
and  finding  the  parameters  of  that  model. .  The  actual  solution  is  in  a 
very  simple  form  and  could  probably  be  implemented  as  an  extension  of 
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state  estimation  techniques.   It  is  by  no  means  the  ideal  answer  to  how 

the  outside  world  should  be  modeled,  but  it  does  allow  a  prediction  of 
how  the  outside  world  will  interact  with  the  systen  of  interest. 
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APPETDEC  A.      M3NLIIEAR  LOAD  FLCW  EQUATIONS 

For  the  pi  transmission  lire  model  shown  in  Figure   2.1,   the   current 
flowing   into  the    line   from  bus   i  is 

^ik  ■  £i  ^ik +  (Ei  -  \'  ^ik  <*•« 

when  the  line  is  connected  between  buses  i  and  k.   vs.,  is  the  shunt 

J    lk 

capacitive   admittance,  y.,    is  the  series   admittance,    and  E.    and  E,    are 
the  voltages    at  buses    i  and  k  respectively.      All  are  complex  quantities 
and 

-  Cik 

ySik   =  J>ik   =  *»—  (A'2) 

"J0ik  1  ,.    ,, 

yik     =^ike  =     R.v+   jX.,  (A'3) 

lk       J    lk 

J'Di 
E        =  E     e     x  (A.  A) 

jDk 
Ek     =  Ek  e     K  (A. 5) 

where  R.v.»  X.,  ,  and  ys.,  are  the  resistance,  reactance,  and  capacitive 
susceptance  of  the  transmission  line  between  buses  i  and  k.  D  and  D, 
are  the  voltage  phase  angles  at  buses  i  and  k  as  measured  with  respect 
to   a  reference.     Equation   (A.l)   can  be   rewritten 

IL1!t  .  E.(^.k  ♦  yik)   -  y.kEk  (A.6) 

By  Kirchhoff's   current   law,   the  sum  of  currents   entering  the   transmission 
lines   connected  to  bus    i  must  equal  the   current  being   injected    into  the 
bus  by  generators    or  loads    (positive   for  generators    and    negative   for 
loads).      I     is   the   current   injected   into  bus    i. 
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\   -    ?  \y   -  *±    2     SJn  *  ylk)    -  S    ylkEk  (A.  7) 


Let  Y,         be   a  complex  bus    admittance  matrix  whose  diagonal  element   in 
row  i   is   the  sum  of  shunt   and  series    admittances   of  all   lines   connected 
to  bus    i 

^    _  -je.. 

(Y,      )..    =    ZL,    (vs.,    +  v..)  =  Y..   =  Y.,  e  (A. 8) 

-bus  n         ,        J    lk       Jiky  n         ii  v        ' 

k 

and  whose  off  diagonal  ik  and  ki  elements  are  the  negative  value  of  the 

series  admittance  of  the  line  between  buses  i  and  k  (or  zero  if  there  is 

no  line  between  buses  i  and  k). 

SbuAk  ■  -yik  =  ?ik  ■  Yik  e"°eik  (A-9> 

Y,    is  a  symmetric  matrix.  Then  for  complex  vectors  L    and  E, 
-bus      J  r  -^bus     -bus 

+  Vi 

whose  i   elements  are  the  current  injected  into  bus  i  and  the  voltage 
at  bus  i  respectively,  equation  (A. 7)  can  be  written  in  matrix  form  as 

•^us   -bus  -bus 

When  Kirchhoff's  current  law  is  applied  to  the  entire  network,  it  is 
obvious  that  for  a  system  with  K+l  buses,  the  current  injections  into 
K  of  the  buses  determine  the  current  injected  at  the  (K+l)   bus. 
Similarly,  the  voltage  angles  are  relative  to  one  another  so  one  bus 
must  be  designated  the  reference  bus  for  voltage  angle  whereas  the 
reference  for  voltage  magnitude  is  ground. 

Real  and  reactive  power  injected  into  bus  i,  P.  +  jQ.,  in  terms  of 
the  current  injected  and  bus  voltage  is 

P.  +  n±   =   E.  I*  (A.  11) 

or  P.  -  jQi  **^\  (A.  12) 


_Q4- 


—   * 


where  E.      is   the   complex  conjugate   of  voltage   at  bus    i.      Using 
equation    (A. 7)    alo 
equation    (A.l?)    is 


equation    (A. 7)    along  with    (A. 8)    8nd    (A. 9)   where  Y..    is  defined, 


P.    -  jQ.    =  E.      y    E,    Y.. 
li  l     A-t      k     lk 


K+l 

z 

k=l 


-jD    ui        jd  -jelk 

2  Ek  e     k  Y.k  e        lk  (A.  13) 


=  v 

k"*l 


K+i  -j(e.,+D-D) 

P.    -   JQ.    =    2    E-EvY-v  e  U.H) 

l       d    i       ,*-i      i  k  lk 
k=l 

The   load   flow  equation,    (A.  14.),   can  be  used   for   networks   consisting  of 
transmission  lines    and   fixed  tap  transformers.      However  for  other  cir- 
cumstances such  as  voltage  controlled  buses,   tap  changing  tr a 21s formers , 
and  phase   shifting  transformers   the   equation  must  be  modified.     Stagg 
and  El-Abiad  1.10  J  is   one  reference  for  these   and  other  variations. 
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APFE^DH  3.      LINEAR  LOAD  FLOW  EQUATIONS 

The   linear   approximation  for  real  power   flow  in  a  transmission 
line  between  bus    i  and  bus  k  is 

where  D.    and  D,    are  the  voltage  phase   angles   at  buses   i   and  k  as 
measured   with  respect  to   a  reference,  E     and  E,    are  the   nominal  values 
of  the  bus  voltage  magnitudes,    and  y.,    is   the  transmission  line   admit- 
tance magnitude.      This   approximation  is   subject  to  the  conditions 

1.  D.    -  D,    is   small. 

l         k 

2.  The  bus  voltage  magnitudes   do   not   vary  much   from  their 
nominal  values. 

3.  The  ratio  of  transmission  line  resistance  to  reactance  is 
small. 

For  convenience   it  will  be   assumed  that  E.    =  E      =  1.0    (or  equivalently 
y.,    is    normalized   for  E.    and  E,  )   so  that  equation   (B.l)    is 

PLik  -  yik(Di  -  V  =Vik-^ikDk  (E-2> 

Assuming  no  losses  within  the  bus,  by  continuity  of  power  the  power 
being  injected  into  a  bus  (positive  for  generators  and  negative  for 
loads)  is  equal  to  the  algebraic  sum  of  the  power  flows  into  the  trans- 
mission lines  conrected  to  that  bus.   Let  P.  be  the  power  injection  at 
bus  i,  then 

pi=  ?  PLik=Di  2.  ylk-  2  yik\  (b.3) 

k  k        k 

Define  P  to  be  a  vector  whose  i   element  is  the  power  injected  into  bus 
i,  and  define  D  to  be  a  vector  of  corresponding  bus  voltage  angles. 
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Then  relation  (B.3)  becomes 


P  =  Y  D 


(3.A) 


Ey  comparison  with  (B.3)  it  can  be  seen  that  the  elements  of  Y  are  such 

that  the  diagonal  element  of  row  i  is  the  sum  of  admittance  magnitudes 

of  the  transmission  lines  which  are  connected  to  bus  i,  Y. .  =  A.  v.,  . 

'        11         ,      Jik 

k 

The   off  diagonal  elements    ik  and  ki  are   negative  value   of  the   admittance 
magnitude   of  the   line  between  bus    i  and  bus   k,     ¥.,    =  Y,  .    =  -v.,  . 
Therefore  the  matrix  Y   is   of  the   form 


Y   = 


^  ylk 

k   1K 

"y12 

"yl3 

-   -   - 

"y12 

^y2k 
k     dK 

"y23 

-  _   - 

"y13 

"y23 

Zy3k 

_   _   _ 

(3.5) 


Because  each  bus  is  usually  only  connected  to  a  few  of  all  the  other 
buses,  Y  is  normally  a  sparce  matrix.   It  is  readily  observable  that 
the  sum  of  all  elements  of  any  row  or  column  of  Y  is  zero.  The  bus 
voltage  phase  angles  must  be  measured  with  respect  to  a  common  reference 
so  one  bus  must  be  designated  as  the  reference  bus  where  voltage  angle 
is  specified.   Also,  because  the  linear  model  Is  lossless,  by  continuity 
of  power,  the  algebraic  sum  of  elements  of  P  is  zero.   Hence  one  row  of 
equation  (3.4.)  is  redundant.   If  bus  i  is  chosen  as  reference  for  vol- 
tage phase  angle,  then  let  D,    be  the  vector  D  with  element  i  deleted, 

let  P,    be  the  vector  P  with  element  i  deleted,  and  let  Y.    be  the 
-bus  -  -bus 


r. 


atrix  Y  with  row  i  and  column  i  deleted.   Then  equation  (B.4)  is 


P.    =  Y,    D. 
-bus   -bus  -bus 


(B.6) 
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P,    and  D,    can  be  separated  into  vectors  PI  8nd  DI  for  Internal 
-bus     -bus         -  —     — 

System  (IS),  F3  and  DB  for  Boundary  System  (3S),  and  PX  and  DX  for 
External  System  (XS).   Similarly  separating  Y,    into  submatrices 


II 

PB 

PX 


YII  YB  Yn 

YB'  YB3  YBX 

YIXT  YBX'  YXX 


DI 
DB 
DX 


(3.7) 


However  by  definition  there  are  no  transmission  lines  between  buses  of 
IS  and  XS  so  YDC  =  0. 


ZI 

III 

YB 

0 

21 

PB 

= 

YI3f 

YBB 

YBX 

D3 

PX 

0 

YBX' 

YXX 

DX 

(B.8) 


The    inverse   of  the  bus   admittance  matrix  Y,         is  the  bus    impedance 

-bus 


matrix  Z,        such  that 
-bus 


D 


-bus       ^bus  -bus 


(B.9) 


or 


DI 

DB 
DX 


ZII 


ZIX1 


ZTB 
2B3 
ZBX' 


ZIX 

PI 

ZBX 

PB 

zxx 

PX 

_       _ 

Z,         in  terrr:s   of  submatrices   of  Y,         can  be    found  by  solving 
-bus  -bus  "  ° 

Z,         Y,         =1  for  Z,         where    I  is  the    icentitv  matrix. 
^us  -ous        —  -bus  — 


ZII         ZIB 


ZIX 


ZJB»        ZBB  ZBX 

ZIX'        ZBX'        ZXX 


YII         YIB  0 

YB'        YBB         YBX 
0  YBX '        YXX 


r 


I 

0 

0 

0 

I 

0 

0 

0 

T 

(B.10) 


(3.11) 


-88- 

Multiplying  the  middle  row  of  Z,         bv  the    left  column  of  Y,         (here 
y  ^bus     -  -ous 

rows   and   columns   will  refer  to  row  and  column  sub-matrices .) 

ZB'  Yrj+   Z3B  YJB'    =  0  (B.12) 

and  solving  for  ZJB' 

ZB'    =  -Z3B  YB'   YII"1  (3.13) 

Multiplying  the  middle   row  of  Z,         by  the  right  column  of  Y,         and 
solving   for   ZBX 

Z33  YBX  -f   ZBX  YXX   =  0  (B.U) 

ZBX   =  -Z33  YBX  YXX"1  (B.15) 

Mult  inlying  the  middle   row  of  Z,         by  the  middle  column  of  Y,         and 
-  -bus  -bus 

solving   for   ZBB  with  the   use   of    (B.I3)    and    (3.15) 

ZB'   YB  +   ZBB  YB3  +   ZBX  YBX'    =   I  (B.16) 

-ZBB  YJB'   YII"1  YB  +   Z3B  YB3  -   ZB3  YBX  YXX^YBX'    =   I  (B.17) 

ZB3    =    (-YB'   YII"1  YB  +  Y33    -  YBX  YXX"1  YBX')"1  (3.18) 

Multiplying  the  top  row  of  Z,        by  the   left  column  of  Y,  and  solving 

<)  LiS  «J  <-lo 

for  ZII  with  the  use   of    (B.I3) 

ZII  YII  +  ZB  YB'    =   I  (B.19) 

ZII  =  YII"1   -  ZB  YB'   YII"1  (3.20) 

ZII  =  YII"1  +  YII"1  YB   ZBB  YJB'   YII*1  (B.21) 

Multiplying  the   too  row  of  Z,        by  the  right  column  of  Y,  and  solving 

^  J  ~QUS  — OUS 

for   ZLX  with  the   use   of    (B.I3)   again 

ZB  YBX  +   ZH  YXX  =  0  (3.22) 


ZLX   =   -  ZIB  Y3X  YXX"1 


ZDC  =  Y 11" -4 13   Z33  YBX  YXX"1 


Let  matrices  YEQI.   YEQX,   AJJI,    and   AQX  be   defined  by 
YSOT  =  -YD'   YIl"1  YIB 


YEQX   =  -YBX  YXX"1  YBX' 


-1 


AQI  =  -YE'   YII 
AQX   =  -YBX  YXX"1 


Then  equations   for  ZB3,    ZII,    ZJ3,    ZJX,    and   ZBX  can  be  written 
Z3B    =    (YEQI  +  YBB  +  YEQX)"1 
ZII   =  YII"1  +   AQI'    ZB3  AQI 


ZB    =  AQI'    Z33 


ZIX   =  AQI'    Z33  AQX 


ZBX  =  ZB3  AQX 


Using    (B.IC),   the   equations   for  DI  and  DB   are 


(B.23) 
(3.24) 

(3.25) 
(3.26) 
(3.27) 
(B.28) 

(B.29) 
(B.30) 
(3.31) 
(3.32) 


DI 


B3 


ZII 

Z33~ 

21 

+ 

ZIX 

ZIB' 

ZBH 

fb 

ZBX 

PX 


(3. 3  U) 


or  using    (3.29)  through    (B.30) 
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DI 


DP 


YII"1 


0 

+ 

^AQI«" 

0 

I 

ZBB  [A^I        l"f 


PI 


PB 


AOI' 


Z3?    AOX  PX 


(3.35) 


From  equation   (3.35)  the  key  role  of  ZBB   is  obvious. 
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AFPE2DH   C.      NONLINEAR  KODEL  IDENTIFICATION 


To  avoid   notation  clutter   in  this    appendix,    all  voltage   8nd   current 
vectors    and   admittance  matrices   will  be  understood  to  be   complex  vectors 
and  matrices   even  though  there   is    no  bar  over  the  symbol    (for  a  complex 
quantity)   or  bar  under  the   symbol    (for   a  vector   or  matrix).      For  example 
E   and    I  are  complex  vectors. 

The   complex  matrix  equation  for  relating  bus  voltages  to  bus   cur- 
rent  injections   is 


dus         bus     bus 


(CI) 


L         and  E,         can  be   divided   into  separate  vectors    for   Internal  System 
(II  and  EI  respectively),    for  Boundary  System    (IB   and  EB),    and   for 
External  System    (IX   and  EX).     Dividing  Y.         into  the   corresponding 
submatrices   and   noting  that  by  definition  of   IS,  BS ,    and  XS   there   is    no 
immediate  coupling  between   IS   and  XS ,   equation   (C.l)   can  be  written 


II 

IB 
IX 


YII         YIB  0 

Y  IB »        YBB  Y3X 

0  YBX '        YXX 


EI 

EB 
EX 


(C2) 


The   only  coupling  between  OS   and  XS   is   due  to  YBB   which  can  be  separated 
into  two  matrices.      That  belonging  to  OS  will  be  designated  YBOS,    and 
that  belonging  to  XS   will  be  designated  YTI,   so  that  YBB    =  YBOS   +  YTL. 
Separating  equation    (C.2) 


-92- 


II 
IB 
EC 


YII         YIB  0 

YT3'        Y30S        0 

0  0  0 


EI 

EB 
EX 


0  0  0 

0         YTL         YBX 
0         YBX '        YXX 


EI 

E3 
EX 


(C3) 


Removing  the   equation  for  OS 


II 


E 


YII         YIB 


v  j3r        YB  <~S 


EI 


E3 


0 


YTL  EB  +  YBX  EX 


(C4) 


The  quantity  YTL  EB  +  YBX  EX  is  a  vector  of  line  current  leaving  the 
buses  for  3S  over  the  tie  lines  between  BS  and  XS.  Designating  this 
the  tie   line  vector   TIL 


II 


IB-LTL 


YII         YTB 


YI3'        Y30S 


EB 


TTL  =  YTL  EB  +  YBX  EX 


(C5) 


(C6) 


Removing  the  botton  row  from  equation    (C.3) 


LX   =  YBX'   E3  +  YXX  EX 


and   solving   for  EX 


EX   =  -YXX"1  YBX'   E3  +  YXX""1   DC 


Substituting  equation    (C.8)    into    (C.6) 


TIL  =  YTL  EB   -  YBX  YXX"1  Y3X'   E3  +  YBX  YXX-1   LX 


(C7) 


(C8) 


(C9) 


Define  an  equivalent  external  bus  admittance  matrix  YEQX  and  an  equi- 
valent external  complex  bus  current  injection  vector  IEQX  by 
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YEQX   =  -YBX  YXX"1  YBX1 


where 


IEQX   =  AQX    IX 


AQX   =  -YBX  YXX 


-1 


Then  equation   (C.6)  becomes 

TTL  =    (YTL  +  YEQX)EB   -   IEQX 

and   if  equation   (C.13)    is   substituted    into  equation   (C.$),   since 
YBB   =  YBCS  +  YTL 


(CIO) 
(Cll) 
(C.12) 

(C.13) 


II 


TBB+IEQI 


YII         YB 


YB'         (YBB+YEQX) 


EI 


EB 


(CH) 


In  the  same  vay,    an  equivalent  complex  bus   admittance  matrix   and   an 
equivalent  complex  bus   current   injection  vector   for    IS   as   seen  by  BS 
can  be   found.      The   top  row  of  equation    (C.14.)   is 


II  =  YII  EI  +  YIB  EB 


Solving   for  EI 


EI  =  -YIl"1  YIB  EB  +  YII"1   II 


and   substituting   into  the  bottom  row  of  equation    (C.14.) 


(C15) 


(C16) 


I     VTT--*- 


-1 


IB  +    3EQI  =  -YIB'   YII"     YIB  EB  +  YIB'   YII"      II  +    (YBB  +  YEQX)E3        (C.17) 

Defining   an  equivalent   internal  complex  bus    admittance  matrix  YEQI,   and 
an  equivalent   internal  complex  bus   current   injection  vector    IEQI  by 


YEQI  =  -YIB'   YII"1  YIB 


IEQI  =  AQI   II 


where 


AQI  =  -YIB'   YII 


-1 


(C.18) 
(C19) 
(C.20) 
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equation   (C.17)  becomes 


IEQI  +    IB  +    IEQX   =    (YEQI  +  YBB  +  YEQX)EB 


(C21) 


Just  as  YBB  was  divider!  into  YTL  and  YBC6 ,  so  also  YBOS  can  be  divided 
into  that  belonging  to  IS  designated  as  YEBS ,  with  Y3CS  =  YBB  +  YBBS. 
Then  equation   (C.21)    is 

IEQI  +    IB  +    IEQX   =    (YBIS  +  YEQl)E3  +  YBBS  E3  +    (YTL  +  YEQX)EB  (C.22) 

Thus    far  equations    (C.$),    (C.I3),    (C.H),    and    (C.21)    are  the 
pertinent  equations    in  terms   of  bus   current   injections.      To  place  them 
in  terms   of  bus   power   injections   let   a  diagonal  matrix  whose   elements 
are  the  complex  conjugate  values   of  bus  voltages  be 


[>•] 


El 

0 

0 

0 

•* 
E2 

0     — 

0 

0 

3t 

(Nste  that   [E"i  is   a  diagonal  matrix  whose   diagonal  elements   are   iden- 
tical to  the   elements   of  the  vector  E.)      Also,    let  PI  and   QI,   PB   and 
QB,   and   PX   and  QX  be  the   real  and  reactive  bus   power   injections   of 
IS,  BS,    and  XS.      Then  because   for   any  bus    i,      P.    -  jQ.    =  E.    I. 

PI  -  jQI  =  [EI*1  II  (C.23) 

PB   -  jQB    =    [EB*]    IB  (G.2A) 

PX   -  jQX   =  [EX*]    IX  (C.25) 

FTL  -   jQTL  =   [EB*]    TTL  (C.?5a) 


or 
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-1 


II  =  [EI  J      (PI  -  jQl) 


-1 


n  =  (ex  ]    (px  -  jqx) 

Multiplying  equation   (C.21)  by   [EB    ] 

[E3*](EQI  +    IB  +    IEQX)    =  [EB*](YEQI  +  YBB  +  YEQX)EB 


(C26) 
(C.27) 


(C.28) 


Define  equivalent  internal  and  external  bus  power  injection  vectors  by 


-1 


PEQI  -  jQEQI  =  [EB'~]    IEQI  =  [eb'  ]  AQl[Ef]      (PI  -  jQl) 


-1 


PEQX   -   jQEQX   =  [EB*]    IEQX   =  [EB"]  AQX[EX*1      (PX  -  jQX) 


Then  (C.23)  becomes 


(C.28) 


(C30) 


(PEQI  -   jQEQI)  +    (FB    -   jQ3)  +    (PEQX   -   jQEQX)    =  [EB*]  (YEQI  +  YBB  +  YEQX)EB 

(C31) 


SUMMARY 

For  the   following  a  bar  over   a  symbol  denotes   a  complex  quantity 
while  matrices    and  vectors    are  denoted  by  underlined  symbols.      Using 
equations    (C.5),    (C.23),    (C.24-)*    and    (C.25a)   the    nonlinear   load   flow 
equation  for  OS   is 


PI  "   321 


(FB    -   JOB)    -    (FT_L  -   ,jQTL) 


[EI*]        0 

o     [IS*] 


YII         YB 


YIB  Y30S 


ET 


EB 


(C32) 


which,    if  YEQX   and   PEQX"  -   jQEQX   are   known,  becomes 
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u  -  m 


(FB    -   JOB)   +    (PEOX    -   .jQEQX) 


[ei*]      0 


0       [E3SJ 


yii 


YIB 


YJ3         (Y5E+YF0X) 


EB 


(C33) 
The    nonlinear   analog   of  the    linear   load   flow  tie   line  power  flow  model 
for   identifying  YE  OX   is 

FTL  -   jOTL  =  [EB*]    (YTL  +  YEQX)E3   -    (FrOX   -   jQEQX)  (C.3-0 

and   the    nonlinear   analog  of  the    linear  boundary  bus    impedance  model  is 


{(PEOl  -  JO^OI)   +    (FB   -   J23)]  +    (PEQX   -  jQEQX) 
=  [EB*][YEQI  +  YEB  +  YEQX]  EB 


(C35) 


Repeating  the   equations    for  current   injections 


II 


IB    -    TIL 


YII         Y13 


Y  B '        YB1 


EB 


(C36) 


II 


IB  +    IEQX 


YII         YIB 


YB'       (YBB+YEOX) 


EI 


EB 


(C37) 


TTL  =    (YTL  +  YE0X)EB    -    IEQX 


(IEQI  +  JE)  +    IEQX   =    (YEQI  +  YB3  +  YE0X)E3 


(C3*) 
(C39) 
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APPENJLX  D.     MAXIMUM  LIKELIHOOD    WE  NT  IF  IC  AT  ION  EQUATIONS 

Given  a  sequence   of  Input  measurements     u(l),   u(2),...,   u(N)» 
a  sequence   of  corresponding   output  measurements      z(l),    z(2),...,    z(N), 
a  set   of  confidence  coefficients     c(l),   c(2),...,   c(N)     which   indicate 
the   relative  confidence   associated  with  each   input/output   set,    and 
given  the   input/output  relation 

z(n)    =  H  u(n)  +  v(n)  (D.l) 

where  v(n)  is  a  noise  vector  having  a  Gaussian  probability  distribution 

with  zero  mean  and  covariance  —7 — >:  R,  it  is  desired  to  find  the  maximum 

c(n)  -• 

likelihood   estimates   of  H  and  R  subject   to  the   following  conditions: 

1.  The   elements   of  H   are   independent  of  one   another. 

(e.g.    in  general  H  /  H') 

2.  E  ^v(n)^     =  0 

3.  E  £v(n)v'(m)j     =  0     for  n  /  m 
A.      E{v(n)v'(n)}     -     ^y  R 

5.     E  £v(n)  u'(n)j   =  0 
The  probability  distribution  function  for  the  Kxl  vector  v(n)   is 

pfc(n)]    =    [(Tirf  c_K(n)lRl]    2     e   2  (D.2) 

The   likelihood   function  for   z(n)    and  u(n)   given  H  and  R  where   z(n)    and 

u(n)    are   Kxl  vectors,    and  K  and  R  are   KxK  matrices    is 

1 
P[z(n),u(n):  H,Rj   '=    [(2m)K  c"K(n)  \R\]"2     .'W^aM^E]  (D.3) 

where 
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j[z(n),u(n):H,R]  =  |[z(n)  -  H  u(n)]  c  (n)R"1[z(n)  -  H  u(n)]'     (T).A%) 

Using  trjAl  to  mean  the  trace  of  matrix  A  and  using  the  identity 

a'a  =  trja  a'V  where  a  is  a  vector  equation  (D.A)  can  also  he  written 

j[z(n),u(n):H,R]  =  |  tr  fc  (n)R_1[z (n)  -  H  u(n)][z(n)  -  H  u(n)V|   (D.5) 

let  the   sets   Z,    U,    and  V  he   such  that      Z   =  fz(l),    z(2),...,    z(NK  , 
U  =  ^u(l),   u(2),...,   u(N)|,      and     V   =  fv(l),   v(2),...,   v(N)"£  . 
Because  v(n)   is   assumed  to  he  uncorrelated  in  time      (Erv(n)v'  (m)r   =  0 
for   n  /  m)    and  because  v(n)    is   Gaussian,   the   joint  probability  distri- 
bution function  for  V   is   the  product   of  the  distribution  functions    for 
v(l),   v(?),...,   v(N). 

p(V)    -  4v(l)]P[v(2)]...    p[v(N)] 

N 
M  I       4     2  v'(n)c(n)R-1v(n) 

-    MW   IE   c-K(n)]"2     e        *=1  (D.6) 

n=l  J 

The   likelihood    function  of  Z  and   U  given  K  and  R   is   then 

P(Z,U:H,R)    =    [(2^)NK\RlN    J£  c"K(n)     2     e  (D.7) 

n=l  J 

where 

J(Z,U:K,R)    =     \    Zs  trrc(n)R"1[z(n')    -  K  u(n)][z(n)    -  H  u(n)]'HD.S) 
n=l         c  ' 

It  is  desired  to  find  the  values  of  H  and  R  which  maximize  the  likelihood 
function  p(Z,U:H,R).   These  values  of  H  and  R  are  the  maximum  likelihood 

A        A 

estimates  H  and  R.   If  the  logarithm  of  the  likelihood  function  is  maxi- 
mized, then  so  will  the  function.  The  log  likelihood  function  is 
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N 


/n[p(Z,U:H,R)]    =  -  \   NK  /n(2rr)    f  \    g   K  fn  c(n) 


n=l 


-  \   N  ^n/Sl    "  J(Z,U:H,R) 


A   new  function  can  be  defined 


(D.9) 


f(2,U:H,R)    =  -2/n[p(Z,U:H,R)]     -    KK  /n(2rr)    f    2    K  ^  c(n) 

n=l 

=   N  ^n(R|   +   2J(Z,U:H,R) 


(D.10) 


f(Z,U:H,R)   only  has   terms  which   involve  H  or  R,  hence,   minimizing 
f  (Z,U:H,R)   is   equivalent  to  maximizing  the   likelihood   function. 

To  perform  the   necessary  mathematical  manipulations   certain 
matrix  gradient   identities   are  used   from  Athans   and  Schweppe  Ell. 
The  matrix  gradient   for  a  scalar  function  of  a  KxK  matrix  X   is   defined 
for  use  here  as 

'  af  (x)  df  (x) 


<?x 


Cf(x)] 


c?X 


11 


PX 


Kl 


<?f(X) 


<?f(X) 


9X 


1M 


<?X\ 


KM 


(D.ll) 


Two  properties   of  trace   of  a  matrix  are 
tr{A^    =     tr  £a»^ 
tr  [A  B^    =     tr  {b  hi 

The    identities   used   are 


(D.12) 


JltT{-1}     =t 


(D.l/J 
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d 


fj  tr  [A  X  B  X'7    =  3  X'    A  +  B'   X'    A'  (D.15) 

-jj  tr  $A  X"1  b\    =  -X"1  B   A  X"1  (D.16) 


-Jf    ^n|X[    =  X"1  (D.17) 

When  f(Z,U,H,R)    is    at   its  minimum  value,    the  matrix  gradients   of 

f  (Z,U:H,R)   with  respect   to  H  and  R  v/ill  be   zero.      Then  the  maximum  like- 

A  A 

lihood   estimates  H   and  R  will  be  the   solutions    for  K  and  R  of 

-~  f  (Z,U:H,R)    =  0  (D.18) 

and  —  f(Z,U:H,R)    =  0  (D.19) 


<?R 


for 


and 


J 


f(Z,U:H,R)    =   N/n|R|    +   2J(Z,U:H,R)  (D.20) 


(Z,U:H,R)    =    J    S   tr  ^c(n)R_1[z(n)   -Ku(n)][z(n)    -  H  u(n)]|(D.2l) 
n=l 


First,   the  gradient  with  respect  to  H   is 

—^  f  (Z,U:H,R)    =     ~g  2J(Z,U:H,R)    =  0 


=  ~/k    S    trfcWR'^Cn)    -  H  u(n)][z(n)    -  H  u(n)V} 
-  n=l        *" 

N      a  C  l  i 

2,    -JjT     tr  Jc(n)R"±z(n)zMn)    -  c(n)R     z(n)u«  (n)H' 


N 

n=l     '  - 


-  c(n)R"  H  u(n)z'(n)  +  c(n)R~  H  u(n)u'(n)H'2  (D.22) 
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Using   equations    (D.12)    and    (D.I3) 


tr  JcCn^R'^Cn)^'  (n)K'£     =  trfc(n)H  u(n)z'  (n)R_1j- 

=  tr  Jc(n)u(n)z'(n)R"1H,l 
tr  £c  (n)^"1^  u(n)z '  (n)^     =  tr  i*c  (n)u(n)z >  (n)^"3^^ 


(D.23) 
(D.24) 


Then  equation   (D.22)    is 
N 


0  =    2     ~JT{     tr  |-2c(n)u(n)z'(n)R"1H  +  c(n)R_:lH  u(n)u' (n)H'"£      (D.25) 


and   using  equations    (D.1/0   and    (D.I5) 

N 

0=2    L-2c(n^(n)z'(n)E"     +   2c(n)u(n)u'(n)H'R     ] 
n=l 


=  o 


r       N  N  -, 

I-    2  c(n)u(n)z'(n)   +    2    c (n)u(n)u« (n)H'   R" 
L    n=l  n=l  J 


CD. 26) 


Then  unless  R        =  0    (which  would  mean  infinite   noise  covariance) 

N  N 

2T    c(n)u(n)z'(n)    =    2  c (n)u(n)u» (n)H» 
n=l  n=l 

N  N 

H    g   c(n)u(n)u'(n)    =    2  c (n)z(n)u' (n) 
n=l  n=l 

and  the  value  of  H  which  maximizes  p(Z,U:H,R)  is 


(D.27) 


(D.2B) 


A 

i-: 


n  1   r  K  1 

Z  c(n)z(n)u«(n)         £   c(n)u(n)u»  (n) 
n=l  J    Ln=l  J 


-1 


(D.29) 


Taking  the  matrix  gradient  of  f(Z,U:H,R)   with  respect   to  R 
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7r  f(Z.HsH,R)    =  0 


-■£ 


j^nlRl  +     y^  2    tr^cWR^di)    -  H  u»]  [z(n)   -  H  u(n)]j 


CD.  30) 

Using  equations    (D.16)    and    (D.17) 
N 


Q  =  I5"1  -  2T  c(n)R"1[z(n)  -  H  u(n)][z(n)  -  H  u(n)}'  R" 


■  "fMs  -  I  2  c(n)[z(n)  -  H  u(n)][z(n)  -  H  adiJj'U"1       (D.31) 

Again  unless  R   =  0,  using  equation  (D.29)  the  value  of  R  which  mini- 
mizes p(Z,U:H,R)  is 

5=^2   c(n)[z(n)    -  H  u(n)][z(n)    -  H  u(n)]  (D.32) 

n=l 


A  A 


Hence  H  and  R  or  equations    (D.29)    and    (D.32)    are   the  maximum   likelihood 
solutions   to  the  problem. 
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APPENDIX  S.   COMPUTER  SIMULATION  PROGRAM 
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